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ABSTRACT 


The  Joint  Munitions  Effectiveness  Manuals  (JMEM)  were  developed  by  the  Joint 
Technical  Coordinating  Group  for  Munitions  Effectiveness  (JTCG/ME)  to  provide  a  set 
of  data  and  methodologies  that  would  permit  a  standardized  comparison  of  weapon 
effectiveness  across  all  service  communities.  In  recent  years,  the  JMEM  are  being 
integrated  into  a  single  software  program  that  allows  users  to  determine  the  effectiveness 
of  weapon  systems  against  a  specified  target  irrespective  of  the  weapon  delivery  mode. 
As  part  of  the  upgrading  effort,  this  thesis  aims  to  develop  a  program,  written  in  Visual 
C++,  to  automate  the  calculation  of  the  Dilution  of  Precision  (DOP)  associated  with  the 
delivery  accuracy  of  GPS  guided  weapon  systems.  The  DOP  values  generated  by  the 
program  were  compared  with  those  generated  by  commercial  DOP  calculation  software 
for  validation.  Relationship  between  the  Vertical  DOP  and  Horizontal  DOP  as  well  as  the 
effect  of  using  outdated  almanac  information  to  calculate  DOP  values  were  studied.  It 
was  found  that  the  loss  of  one  visible  satellite  could  cause  the  DOP  to  increase  by  as 
much  as  38%. 
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I.  INTRODUCTION 


A.  GLOBAL  POSITIONING  SYSTEM  (GPS) 

1.  Overview  of  the  NAVSTAR  Global  Positioning  System 

The  NAVSTAR  Global  Positioning  System  (GPS)  is  a  satellite-based  navigation 
and  positioning  system  made  up  of  a  constellation  of  between  24  to  32  satellites  or  space 
vehicles  (SV).  The  GPS  was  developed  by  the  United  States  Department  of  Defense  and 
although  the  GPS  was  originally  intended  for  military  applications,  the  United  States 
government  made  the  system  freely  available  for  civilian  use  in  the  1980s.  GPS  works  in 
any  weather  conditions,  anytime  and  anywhere  in  the  world. 

The  system  consists  of  three  segments:  the  space  segment,  the  control  segment, 
and  the  user  segment.  The  United  States  Air  Force  maintains  and  operates  the  space  and 
control  segments. 

a.  Space  Segment 

The  space  segment  consists  of  the  orbiting  satellites.  The  GPS 
constellation  has  a  minimum  of  24  satellites  traveling  on  six  medium  Earth  orbits 
(altitude  about  20,200  km)  of  approximately  55°  inclination  (tilt  relative  to  Earth's 
equator)  and  are  separated  by  60°  right  ascension  of  the  ascending  node  (angle  along  the 
equator  from  a  reference  point  to  the  orbit's  intersection).  Each  satellite  completes  one 
orbit  in  slightly  less  than  12  hours. 

Each  satellite  transmits  its  own  unique  microwave  signals  on  two  different 
L-band  frequencies  that  give  infonnation  on  the  precise  orbit  for  the  satellite  sending  the 
message  (the  ephemeris);  the  approximate  orbits  and  general  health  of  all  satellites  (the 
almanac);  as  well  as  an  ionospheric  delay  model.  All  satellites  broadcast  at  the  same  two 
frequencies,  1.57542  GHz  (LI  signal)  and  1.2276  GHz  (L2  signal).  The  receiver  can 
distinguish  the  signals  from  different  satellites  because  they  are  encoded  with  a  pseudo¬ 
random  number  (PRN)  sequence  that  is  different  for  each  satellite.  The  receiver  knows 
the  PRN  codes  for  each  satellite  and  uses  this  to  reconstruct  the  navigation  message. 
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b.  Control  Segment 

The  control  segment  consists  of  five  monitor  stations  around  the  Earth  that 
maintain  the  satellites  in  their  proper  orbits  through  occasional  maneuvers,  and  adjust  the 
satellite  clocks.  It  tracks  the  satellites,  uploads  updated  navigational  data,  and  maintains 
health  and  status  of  the  satellite  constellation. 

c.  User  Segment 

The  user  segment  consists  of  the  GPS  receiver  equipment,  which  receives 
the  signals  from  the  satellites  and  uses  the  transmitted  information  to  calculate  the 
receivers’  positions,  velocities  and  headings  based  on  the  positions  of  the  satellites  and 
the  time  the  signal  was  transmitted  and  received. 

2.  Workings  of  GPS 

Infonnation  on  the  positions  of  the  satellites  is  transmitted  by  the  satellites  and 
can  be  calculated  relative  to  a  set  of  coordinates  that  are  Earth  centered,  Earth  fixed 
(ECEF)  (see  Figure  1). 


Z 


Figure  1.  Earth  Centered,  Earth  Fixed  (ECEF)  Frame  [After  1] 
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The  range,  R ,  between  the  GPS  receiver  and  each  satellite  is  measured  by  timing 
the  delay  between  the  transmission  time  of  the  signal  from  the  satellite  and  the 
arrivaltime  of  the  same  signal  to  the  receiver.  At  i.e.,  the  signal  travel  time.  Since  these 
signals  travel  at  the  speed  of  light,  c,  the  range  between  the  GPS  receiver  and  the  satellite 
is  given  by 

R  =  cxAt  (1.1) 

Three  coordinates  in  the  ECEF  frame  define  the  GPS  receiver's  position.  This 
implies  that  three  different  range  measurements  would  suffice  to  detennine  the  position 
of  the  receiver.  However,  clocks  used  in  GPS  receivers  are  not  as  accurate  as  the  atomic 
clocks  in  the  satellites.  As  an  error  of  a  nanosecond  in  time  measurement  would  result  in 
an  error  of  about  0.3  m  in  range  calculation,  each  range  measurement  needs  to  be 
corrected  to  account  for  the  receiver  clock's  inaccuracies.  As  these  measured  ranges  are 
distorted  by  the  relatively  inaccurate  time  keeping  of  the  receiver's  clock,  they  are  known 
as  pseudoranges. 

Since  the  GPS  receiver  clock  error  is  an  unknown  variable  in  addition  to  the  three 
position  coordinates  of  the  GPS  receiver,  a  minimum  of  four  range  measurements, 
instead  of  three,  are  required  to  resolve  the  GPS  time  and  determine  the  position  of  the 
receiver  using  trilateration. 

Trilateration  can  be  described  as  using  the  pseudoranges  to  form  spheres  around 
the  respective  satellites  such  that  the  position  of  the  receiver  lies  within  the  overlapping 
region  of  the  spheres.  It  should  be  noted  that  trilateration  is  different  from  triangulation  in 
that  trilateration  uses  ranges,  while  triangulation  uses  angles  to  detennine  the  position  of 
a  point. 


3.  Sources  of  Errors  in  GPS 

The  above  description  of  how  GPS  works  assumed  that  there  are  no  other  sources 
of  enor  other  than  the  GPS  receiver  clock.  In  reality,  there  are  many  sources,  which  can 
introduce  error  into  the  calculation  of  the  GPS  receiver  position.  These  enors  need  to  be 
accounted  for  in  order  to  mitigate  their  effects  on  the  accuracy.  The  sources  of  errors 
include  the  following: 
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a. 


Atmospheric  Effect 


As  the  satellite  signal  passes  through  the  atmosphere,  its  speed  is  reduced 
as  air  has  a  slightly  higher  index  of  refraction  (about  1 .0003  [2])  causing  the  signal  speed 
to  decrease  by  an  average  of  about  0.03%.  However,  inconsistencies  in  the  atmosphere, 
especially  the  ionosphere,  cause  the  signal  to  slow  in  a  non-uniform  manner.  This  effect 
is  least  when  the  satellite  is  directly  overhead  and  become  greater  for  satellites  near  the 
horizon  since  the  signal  path  through  the  atmosphere  is  longer. 

In  order  to  mitigate  this  effect,  after  the  receiver's  approximate  position  is 
known  using  the  pseudoranges,  a  built-in  mathematical  model  can  be  used  to  estimate  the 
average  amount  of  delay  to  compensate  for  this  type  of  error.  However,  the  model  may 
not  be  able  to  predict  the  full  effects  of  the  ionospheric  delay.  A  more  accurate  way  to 
compensate  for  this  error  is  to  use  both  frequencies  to  measure  the  time  delay. 
Ionospheric  delay  affects  the  speed  of  microwave  signals  differently  depending  on  their 
frequency,  this  is  a  characteristic  known  as  dispersion.  Delays  measured  on  two 
frequency  bands  can  be  used  to  measure  dispersion,  and  this  measurement  can  then  be 
used  to  estimate  the  delay  at  each  frequency. 

Another  way  to  compensate  for  the  ionospheric  error  is  to  compare  the 
GPS-measured  position  with  a  known  surveyed  position.  This  takes  advantage  of  the  fact 
that  the  effects  of  the  ionosphere  generally  change  slowly  and  can  be  averaged  over  time; 
hence,  the  correction  on  the  ionospheric  error  can  be  applied  to  other  GPS  receivers  in  the 
same  general  region.  Satellite  Based  Augmentation  Systems  (SBAS)  such  as  WAAS 
(available  in  North  America  and  Hawaii),  EGNOS  (Europe  and  Asia)  or  MSAS  (Japan) 
transmits  the  ionospheric  correction  data  via  satellite,  while  Ground  Based  Augmentation 
Systems  (GBAS)  transmits  the  correction  data  via  ground  radio  transmitter  directly  to  the 
GPS  receiver. 

Humidity  in  the  troposphere  also  results  in  errors  similar  to  ionospheric 
delay.  However,  this  effect  is  more  localized,  changes  more  quickly  than  ionospheric 
effects,  and  is  not  frequency  dependent.  These  characteristics  make  tropospheric  effects 
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more  difficult  to  measure  and  compensate  compared  to  ionospheric  effects.  Typically, 
error  in  pseudorange  caused  by  ionospheric  effects  are  about  ±  5  m  while  the  troposheric 
effect  is  about  ±  0.5  m  [3]. 

b.  Multipath 

When  the  GPS  signal  is  reflected  off  objects  such  as  tall  buildings,  the 
travel  time  of  the  signal  before  it  reaches  the  receiver  increases.  This  results  in  multipath 
errors.  Various  techniques  have  been  developed  to  mitigate  multipath  errors.  For  long 
delay  multipath,  the  receiver  itself  can  be  programmed  to  recognize  the  wayward  signal 
and  discard  it.  For  shorter  delay  multipath  from  the  signal  reflected  off  the  ground, 
specialized  antennas  may  be  used  to  reduce  the  power  received  by  the  antenna  from  such 
reflected  signals.  Short  delay  reflections  are  harder  to  distinguish  from  routine 
fluctuations  in  atmospheric  delay. 

Multipath  effects  are  less  severe  in  moving  vehicles  as  solutions  using 
reflected  signals  quickly  fail  to  converge  and  only  the  direct  signals  result  in  stable 
solutions.  Typically,  error  in  pseudorange  caused  by  multipath  effect  is  about  ±  1  m  [3]. 

c.  Ephemeris  and  Clock 

The  satellites  transmit  ephemeris  data  (data  on  their  precise  orbits)  every 
30  seconds,  but  the  data  itself  may  be  up  to  two  hours  old.  Although  data  up  to  four  hours 
old  is  considered  valid,  it  may  not  indicate  the  satellite's  actual  position. 

The  satellite's  atomic  clocks  encounter  noise  and  clock  drift  errors.  While 
the  navigation  message  contains  corrections  for  these  errors  and  estimates  of  the  accuracy 
of  the  atomic  clock,  they  are  based  on  observations  done  at  the  monitor  stations  and  may 
not  indicate  the  clock's  actual  state.  However,  these  errors  are  typically  small.  Error  in 
pseudorange  caused  by  ephemeris  error  is  about  ±  2.5  m  while  clock  error  is  about  ±  2  m 
[3]. 
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d.  Relativity 

According  to  the  theory  of  relativity,  the  clocks  on  the  satellites  are 
affected  by  their  speed  (special  relativity)  as  well  as  their  gravitational  potential  (general 
relativity).  Due  to  the  weaker  gravitational  field  at  the  GPS  orbit,  general  relativity 
predicts  that  time  speeds  up  by  about  45.9  ps  per  day.  On  the  other  hand,  special 
relativity  predicts  that  time  slows  by  about  7.2  ps  per  day  due  to  the  orbital  speed  of  the 
satellite.  Hence,  the  total  effect  is  that  time  on  the  satellite  speeds  up  by  about  39  ps  per 
day. 

Since  accurate  time  keeping  is  central  to  the  accuracy  of  GPS,  this 
discrepancy  has  to  be  accounted  for  and  this  is  done  by  giving  the  frequency  standard  on 
board  each  satellite  a  rate  offset  prior  to  launch,  making  it  run  slightly  slower  than  the 
desired  frequency  on  Earth;  specifically,  at  10.22999999543  MHz  instead  of  10.23  MHz. 

e.  Sagnac  Distortion 

Sagnac  distortion  is  caused  because  GPS  time  is  defined  in  an  inertial 
frame  while  observations  are  processed  in  an  ECEF  frame.  A  Lorentz  transformation  is 
applied  to  convert  from  the  inertial  frame  to  the  ECEF  frame  and  the  resulting  correction 
on  the  signal  travel  time  has  opposite  algebraic  signs  for  satellites  in  the  Eastern  and 
Western  celestial  hemispheres.  Although  the  effect  is  small,  neglecting  it  will  produce  an 
east- west  error  of  about  a  few  hundreds  of  nanoseconds,  or  tens  of  meters  in  position  [3], 

f  Selective  Availability 

Selective  Availability  (SA)  is  a  feature  in  GPS  that,  when  enabled,  can 
introduce  intentional  random  errors  of  up  to  a  hundred  meters  into  the  civilian  navigation 
signals  with  the  intention  of  limiting  accurate  positioning  capability  to  the  United  States 
military  and  other  authorized  users. 

However,  this  feature  was  turned  off  in  May  2000  following  an  executive 
order  from  the  United  States  President  Bill  Clinton  to  set  the  SA  error  to  zero  by  2006. 
This  allowed  civilian  applications,  such  as  the  aviation  industry,  to  take  advantage  of  the 
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highly  accurate  navigation  signals.  In  Sep  2007,  the  United  States  Department  of  Defense 
announced  that  future  GPS  satellites  would  no  longer  support  SA,  thereby  making  the 
policy  permanent. 

g.  Jamming 

Having  travelled  about  20,200  km  from  the  satellites  to  Earth,  GPS  signals 
received  by  GPS  receivers  on  Earth  tend  to  be  relatively  weak.  Hence,  it  is  easy  for  other 
sources  of  electromagnetic  (EM)  radiation  to  overpower  the  GPS  signals,  making 
acquiring  and  tracking  the  satellite  signals  difficult  or  impossible.  These  sources  of  EM 
radiation  can  occur  naturally  or  made  artificially. 

An  example  of  naturally  occurring  EM  radiation  capable  of  disrupting 
GPS  reception  is  solar  flare.  Solar  flares  are  explosions  in  the  Sun  that  produces  strong 
EM  radiation  that  has  the  potential  to  disruption  satellite  communication.  Other  examples 
include  naturally  occurring  geomagnetic  storms,  found  mainly  near  the  poles  of  the 
Earth's  magnetic  field  as  well  as  interference  from  Van  Allen  Belt  radiation  when  the 
satellites  pass  through  the  South  Atlantic  Anomaly. 

An  example  of  artificial  source  is  man-made  jammers,  which  typically 
emit  strong  EM  radiation  to  overpower  the  actual  GPS  signal.  These  signals  can  interfere 
with  GPS  receivers  when  they  are  within  radio  range  or  line  of  sight. 

h.  Number  of  Visible  Satellites 

The  GPS  constellation  is  designed  to  have  at  least  six  satellites  above  any 
part  of  the  Earth  at  any  one  time.  However,  reception  of  the  signals  from  these  satellites 
may  be  blocked  by  nearby  obstacles  such  as  buildings,  terrain  and  dense  vegetation 
making  position  calculations  less  accurate.  In  the  worst  case,  all  signals  may  be  blocked 
making  position  calculation  impossible.  In  general,  the  higher  the  number  of  visible 
satellites,  the  better  the  accuracy. 
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4. 


Dilution  of  Precision 


The  dilution  of  precision  (DOP)  also  contributes  to  the  accuracy  of  the  GPS 
calculations  but  not  in  a  direct  manner.  Mathematically,  DOP  is  the  ratio  between  the 
standard  deviations  of  a  specified  parameter  and  the  pseudorange.  For  example,  Vertical 
DOP  is  the  ratio  between  the  standard  deviation  of  the  vertical  component  (altitude)  of 
the  GPS  receiver  and  the  standard  deviation  of  the  pseudorange.  For  parameters  that 
involve  more  than  one  variable  such  as  Geometric  DOP,  the  ratio  is  between  the  root  sum 
square  of  the  standard  deviation  of  the  variables  (x,  y,  z  coordinates  and  time)  and  the 
standard  deviation  of  the  pseudorange. 

Physically,  DOP  describes  the  geometric  strength  of  the  visible  satellites' 
configuration  on  the  GPS  accuracy.  Ideally,  the  visible  satellites  should  be  located  at 
wide  angles  relative  to  each  other.  The  geometry  of  such  satellite  configuration  is  said  to 
be  strong  and  the  DOP  values  are  low.  Conversely,  if  the  visible  satellites  have  small 
angular  separation,  the  satellites'  configuration  has  weak  geometry  and  the  DOP  values 
are  high. 

Figure  2  shows  a  scenario  where  a  GPS  receiver  measures  the  pseudoranges  of 
two  satellites.  Keeping  the  error  of  the  range  measurement  constant  in  both  cases,  the 
case  on  the  left,  with  larger  angular  separation  between  the  two  satellites,  shows  that  the 
area  of  uncertainty  on  the  position  of  the  receiver  is  smaller  than  the  case  on  the  right. 
The  DOP  value,  hence,  may  be  understood  as  the  'dilution'  factor  on  the  accuracy  of  the 
original  measurement. 
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Pseudorange 


Low  DOP 


High  DOP 


Figure  2.  Effect  of  Satellite  Geometry  on  Dilution  of  Precision  [After  13] 
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Table  1.  Meaning  of  DOP  Values  [From  4] 


DOP  Value 

Rating 

Description 

1 

Ideal 

This  is  the  highest  possible  confidence  level  to  be  used  for 

applications  demanding  the  highest  possible  precision  at  all 

times. 

2-3 

Excellent 

At  this  confidence  level,  positional  measurements  are 

considered  accurate  enough  to  meet  all  but  the  most  sensitive 

applications. 

4-6 

Good 

Represents  a  level  that  marks  the  minimum  appropriate  for 

making  business  decisions.  Positional  measurements  could  be 

used  to  make  reliable  in-route  navigation  suggestions  to  the 

user. 

7-8 

Moderate 

Positional  measurements  could  be  used  for  calculations,  but  the 

fix  quality  could  still  be  improved.  A  more  open  view  of  the 

sky  is  recommended. 

9-20 

Fair 

Represents  a  low  confidence  level.  Positional  measurements 

should  be  discarded  or  used  only  to  indicate  a  very  rough 

estimate  of  the  current  position. 

21-50 

Poor 

At  this  level,  measurements  are  inaccurate  by  as  much  as  300 

meters  with  a  6-meter  accurate  device  and  should  be  discarded. 

Table  1  gives  a  description  and  the  meaning  for  various  DOP  values.  It  should  be 
noted  that  DOP  values  of  less  than  1  are  possible  i.e.,  the  accuracy  of  the  calculated 
position  (via  trilateration)  could  be  higher  than  the  accuracy  of  the  pseudoranges 
measured  by  the  GPS  receiver.  When  there  are  sufficient  visible  satellites  with  wide 
angular  separation,  it  is  possible  for  the  region  of  uncertainty  on  the  position  of  the 
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receiver  to  be  reduced  to  the  point  that  it  is  smaller  than  the  uncertainty  of  the  individual 
pseudorange  measurement.  This  is  inherent  in  the  trilateration  method. 

B.  JOINT  MUNITIONS  EFFECTIVENESS  MANUALS  (JMEM) 

1.  Introduction  to  JMEM 

The  Joint  Technical  Coordinating  Group  for  Munitions  Effectiveness  (JTCG/ME) 
produced  the  Joint  Munitions  Effectiveness  Manuals  (JMEM)  in  order  to  provide  a  set  of 
data  and  methodologies  that  standardize  weapon  effectiveness  calculations,  thereby 
facilitating  comparisons  of  weapon  effectiveness  between  different  communities  in  the 
military. 

JMEM  include  detailed  information  on  the  physical  characteristics  and 
performance  of  weapons  and  weapon  systems;  descriptions  of  the  mathematical 
methodologies  that  employ  these  data  to  generate  effectiveness  estimates;  software  that 
permit  users  to  calculate  effectiveness  estimates;  and  pre-calculated  weapon  effectiveness 
estimates.  They  are  used  by  all  services  in  United  States  as  well  as  NATO  and  other  allies 
to  plan  operational  missions,  support  training  and  tactics  development,  and  support  force- 
level  analyses. 

In  the  past,  JMEM  were  volumes  of  orange  covered  manuals  but  over  time, 
computer  programs  were  used  to  supplement  the  manuals  to  allow  faster  computation  as 
well  as  more  realistic  (but  also  more  complex  and  computationally  intensive)  models  to 
be  used.  As  computers  became  more  affordable  and  widely  used,  the  paper  version  of 
JMEM  gave  way  to  CD  versions.  Beginning  in  2007,  all  JMEM  weapon  effectiveness 
products  are  integrated  into  a  single  program  called  the  JMEM  Weaponeering  System 
(JWS).  This  is  a  target-oriented  program,  which  allows  users  to  determine  the 
effectiveness  of  weapon  systems  against  a  specified  target  regardless  of  the  weapon 
delivery  mode. 

2.  Use  of  JMEM  in  Mission  Planning 

During  the  planning  of  offensive  missions,  it  is  important  for  military  planners  to 
know  about  the  delivery  accuracy  as  well  as  the  effectiveness  of  selected  weapon  system 
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against  a  specified  target.  This  will  allow  the  planners  to  estimate  the  number  of 
weapon  systems  required  to  destroy  the  target. 

The  effectiveness  of  a  weapon  system  against  a  specified  target  is  quantified  by 
the  effectiveness  index  or  lethal  area.  This  is  dependent  on  factors  such  as  the  defined  kill 
criterion;  the  physical  and  geometrical  configuration  of  the  target  and  its  critical 
components;  nature  of  the  weapon  system;  and,  the  damage  required  on  the  critical 
component(s)  to  achieve  the  desired  kill  criterion. 

The  delivery  accuracy  of  a  weapon  system,  on  the  other  hand,  is  quantified  by  the 
distribution  of  the  weapon  system’s  impact  points  such  as  deflection  error  probable  (DEP) 
and  range  error  probable  (REP).  The  weapon  system’s  accuracy  is  dependent  on  many 
factors  such  as  target  acquisition  error,  tracking  error,  etc. 

3.  Effect  of  DOP  on  Mission  Planning 

For  the  case  of  GPS-guided  weapon  systems,  the  error  associated  with  GPS  in 
detennining  the  position  of  the  target  and  the  weapon  is  one  of  the  most  important 
sources  of  errors  to  take  into  account  when  calculating  the  delivery  accuracy  of  the 
weapon  system.  As  DOP  has  the  effect  of  amplifying  the  original  GPS  measurement 
errors,  any  mission  planners  who  intend  to  use  GPS-guided  weapon  systems  would  need 
to  know  the  DOP. 

Since  DOP  is  derived  from  the  configuration  of  the  visible  satellites,  it  varies 
depending  on  the  time  and  position  of  the  target.  For  mission  planners,  that  variation  in 
DOP  could  make  a  difference  to  the  effectiveness  of  an  offensive  mission  and  the 
potential  collateral  damage. 

C.  MOTIVATION  OF  THESIS 

As  part  of  the  effort  to  integrate  the  JMEM  into  a  single  software  program,  this 
thesis  aims  to  develop  a  program  to  automate  the  calculation  of  DOP  associated  with  the 
delivery  accuracy  of  GPS-guided  weapon  systems. 
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D.  OBJECTIVES  OF  THESIS 

This  thesis  aims  to  develop  a  program  in  Visual  C++  to  automate  the  calculation 
of  DOP  using  the  almanac  file  in  SEM  format  (,al3).  In  addition,  the  thesis  also  aims  to 
study  the  effect  of  the  using  outdated  almanac  data  in  the  DOP  calculation  and  find  out 
the  relationship  between  the  Horizontal  DOP  and  Vertical  DOP. 
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II.  THEORY 


A.  POSITIONS  OF  GPS  SATELLITES 

The  position  of  satellite  at  any  given  instant  can  be  calculated  from  the  ephemeris 
of  the  satellite.  The  almanac  is  a  practical  and  convenient  source  to  get  the  ephemeris  of 
all  the  satellites  in  the  constellation.  Although  the  almanac  only  gives  the  rough 
ephemeris,  the  accuracy  is  good  enough  for  calculating  DOP  values.  Moreover,  the 
almanacs  are  posted  ahead  of  time,  thereby  allowing  planning  to  be  done. 

1.  Almanac  Data 

The  almanacs  are  available  to  the  public  from  the  United  States  Coast  Guard 
website,  http://www.navcen.uscg.gov/GPS/almanacs.htm,  in  the  form  of  a  file  that  is 
updated  almost  daily.  The  almanac  file  is  available  in  two  fonnats,  namely  SEM  (,al3) 
and  YUMA  (.aim).  Examples  of  an  almanac  file  in  SEM  and  YUMA  format  are  shown  in 
Appendix  A. 

The  YUMA  format  is  more  reader-friendly,  but  the  SEM  format  is  more  compact. 
Hence,  in  view  of  efficiency,  the  SEM  format  was  chosen  as  the  input  file  to  calculate  the 
positions  of  all  the  satellites  at  any  given  time.  Useful  information  contained  in  each 
SEM  almanac  file  is  shown  in 
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Table  2.  It  should  be  noted  that  angles  are  given  in  the  tenns  of  number  of  semicircles, 
hence  it  is  important  to  convert  them  to  radians  before  proceeding  with  further 
calculations. 
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Table  2.  Useful  Information  in  Almanac  File  (SEM  Format)  [After  5] 


Information 

Unit 

Description 

Number  of  Records 

records 

The  number  of  satellite  almanac  records 
contained  in  the  file 

GPS  Week  Number,  WN 

weeks 

The  almanac  reference  week  for  all 
almanacs  in  the  file  as  per  ICD-GPS-200 

GPS  Time  of  Applicability, 
TOA 

sec 

The  almanac  reference  time  for  all 
almanacs  in  the  file  as  per  ICD-GPS-200 

PRN  Number 

none 

The  satellite  PRN  number  as  per  ICD- 
GPS-200.  Used  to  identify  individual 
satellite 

Eccentricity,  e 

unitless 

The  satellite  almanac  orbital  "eccentricity" 
as  defined  in  ICD-GPS-200 

Inclination  Offset,  8ik 

semicircles 

The  satellite  almanac  orbital  "inclination 
angle  offset"  as  defined  in  ICD-GPS-200 

• 

Rate  of  Right  Ascension,  Q 

semicircles/ 

sec 

The  satellite  almanac  orbital  "rate  of  right 
ascension"  as  defined  in  ICD-GPS-200 

Square  root  of  Semi-Major 
Axis,  \[a 

m^ 

The  satellite  almanac  orbital  "square  root 
of  the  semi-major  axis"  as  defined  in  ICD- 
GPS-200 

Longitude  of  Orbital  Plane, 
Q0 

semicircles 

The  satellite  almanac  orbital  "geographic 
longitude  of  the  orbital  plane  at  the  weekly 
epoch"  as  defined  in  ICD-GPS-200 

Argument  of  Perigee,  co 

semicircles 

The  satellite  almanac  orbital  "argument  of 
perigee"  as  defined  in  ICD-GPS-200 

Mean  Anomaly  at  Reference 
Time,  M0 

semicircles 

The  satellite  almanac  orbital  "mean 
anomaly"  as  defined  in  ICD-GPS-200 

Satellite  Health 

none 

The  satellite  health  code  expressed  in 
integer  fonn 
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2. 


Calculations 


a.  Gregorian  to  Julian  Date  Conversion 

Before  information  extracted  from  the  almanacs  can  be  used  to  calculate 
the  positions  of  the  satellites,  the  date  and  time  needs  to  be  specified  and  converted  to  the 
same  form  as  those  in  the  almanac,  namely  the  GPS  week  number  (the  number  of  weeks 
since  Jan  6,  1980-  the  reference  start  date  for  GPS)  and  number  of  seconds  of  that  week. 
This  is  done  by  subtracting  the  specified  date  with  the  GPS  reference  start  date.  However, 
this  is  difficult  to  do  using  the  Gregorian  calendar  (the  internationally  accepted  calendar 
used  today,  for  example,  Mar  27,  2009).  Hence,  the  dates  are  processed  in  Julian  dates 
instead.  The  conversion  from  Gregorian  to  Julian  date,  JD  is  shown  below  [6] 

Y  +  M  ~9 

rn  i  jr  r\  jL  I 

JD  =  367  •  7  -  floor(—  ( Y  +  floor (  ^  ))  -  fl oor(- (floor( - ^ — )  + 1) 

„  ,275  -M  „  __  hr  min  sec 

+  floor  ( - )  +  D  + 172 1028.5  + —  + - + - 

9  24  1440  86400 

timezone  +  daylightsaving  leap  sec 
24  +  86400 

In  the  above  equation,  "Y",  "M",  "D",  "hr",  "min"  and  "sec"  represents  the 
year,  month,  day,  hours,  minutes  and  seconds  of  the  specified  local  date  and  time 
respectively;  "timezone"  represents  the  number  of  hours  offset  from  GMT  or  Zulu  Time, 
for  example,  the  offset  for  Eastern  Standard  Time  (North  America)  is  -5  hours; 
"daylightsaving"  represents  an  additional  hour  when  daylight  saving  is  in  effect;  "leap 
sec"  represents  the  number  of  leap  seconds  added  since  Jan  6,  1980.  The  "floor"  operator 
returns  the  integer  value. 

Using  the  above  conversion,  the  Julian  date  for  Jan  6,  1980  is  2,444,244.5 
days.  The  number  of  weeks  since  Jan  6,  1980,  NumWeek,  is  therefore 

JD  -  2444244  5 

NumWeek  =  floorfl- - - - -)  (1.3) 
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If  NumWeek  is  more  than  1024,  minus  NumWeek  by  1024  until  it  is  less 
than  1024.  This  is  due  to  the  GPS  week  rollover  issue  [7]. 

Finally,  the  number  of  seconds  of  that  NumWeek,  NumSec,  is  given  by 

NumSec  =  (JD  -  2444244.5  -  7  •  NumWeek )  •  86400  (1.4) 


b. 


Satellite  Position  in  ECEF  Frame 


With  the  time  specified,  the  positions  of  the  satellites  can  then  be 
calculated  with  reference  to  the  ECEF  frame.  There  are  two  constants  required  in  the 
calculation,  namely  the  WGS  84  value  of  the  Earth's  Universal  Gravitational  Parameter,  p 

=  3.986005  x  10 14  m3/sec2  and  the  WGS  84  value  of  the  Earth's  Rotation  Rate,  = 
7.2921151467  x  10-5  rad/sec.  Using  data  from  the  almanac  file  for  each  satellite,  the 
following  parameters  are  calculated  sequentially  to  obtain  the  satellite  position  in  the 
ECEF  frame  [8]. 

The  Computed  Mean  Motion,  n0  is  given  by 


(1.5) 


The  Time  since  TOA,  tk  is  given  by 

tk  =  (Nam  Week  -  WN)  ■  604800  +  ( NumSec  -  TOA ) 

The  Mean  Anomaly,  Ma  is  given  by 


M,  =  M0  +  n0tk 


(1.6) 


(1.7) 


The  Kelper's  Equation  for  Eccentric  Anomaly  is  shown  below  where  Ek 
needs  to  be  solved  by  iteration  since  the  expression  for  Ek  is  not  explicit. 


Ek  =  Mk  +  e  ■  sin  Ek 


(1.8) 


The  True  Anomaly,  ok  is  calculated  using  the  value  of  Ek  obtained  from 
Equation  2.7  as  follows. 
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(1.9) 


vk  =  tan 


f (Vl  ~e2  smEk) j - e ■  cos Ek) 
1  (cos  Ek  —  e)/(l  —  e  ■  cos  Ek ) 


} 


The  Eccentric  Anomaly,  Ek  is  then  recalculated  using  ok  obtained  from 
Equation  2.8  and  this  new  value  is  used  for  subsequent  equations. 


e  +  cos  o. 


Ek=  cos"(  *  } 

l  +  e-  cos  v,r 


(1.10) 


The  Corrected  Argument  of  Latitude,  uk  is  given  by 


LI  —  Ok  +  CO 


(1.11) 


The  Corrected  Radius,  rk  is  given  by 

rk  =  A(\-e-cosEk) 

The  Corrected  Inclination,  ik  is  given  by 


4  4  sh 


(1.12) 


(1.13) 


The  Satellite  Position  in  Orbital  Plane  is  given  by 


xk  =  rk  cos  uk 


(1.14) 


yk=rksmuk 


(1.15) 


The  Corrected  Longitude  of  Ascending  Node,  is  given  by 


•  •  • 

Q4=Q0+(Q-Q  e)tk-Ele(TOA)  (1-16) 

Finally,  the  Satellite  Position  in  ECEF  Frame  is  given  by 

xk  =  xk  ■  cos  Clk  -  yk  ■  cos  yk  sin  Elk  (1-17) 

yk  =  xk  •  sin  Qk  +  y'k  ■  cos  yk  cos  Qk  (1-18) 

zk =  y'k  ■ sin  4  (!-i9) 
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c.  ECEF  to  ENU  Conversion 

In  order  to  obtain  the  Horizontal  DOP  value  with  reference  to  the  Earth's 
surface,  the  positions  of  the  satellites  need  to  be  converted  to  the  East-North-Up  (ENU) 
coordinates  relative  to  a  local  reference  point  specified  on  the  Earth  surface.  This  position 
corresponds  to  the  latitude  ( (f> )  and  longitude  ( X )  of  the  GPS  receiver's  position. 

There  are  two  constants  required  in  the  calculation,  namely  the  WGS  84 
value  of  the  Earth's  Semi-Major  Axis,  a  =  6,378,137  m  and  the  WGS  84  value  of  the 
Earth's  First  Eccentricity,  el  =  8.181919084266  x  10"2.  Furthermore,  the  local  reference 
point  is  expressed  in  ECEF  coordinates  to  simplify  the  ECEF  to  ENU  conversion  for  the 
satellites'  positions.  The  conversion  from  geodetic  (X,(j),  and  altitude  (alt))  to  ECEF 
coordinates  (x,  y,  z)  is  shown  below  [9],  It  should  be  noted  that  for  the  local  reference 
point,  alt  is  zero  since  it  is  specified  to  be  on  the  Earth's  surface. 

The  Prime  Vertical  Radius  of  Curvature,  N  is  given  by 

N=  U  (1.20) 

^l-e2  sin2  <j) 

The  conversion  from  Geodetic  to  ECEF  Coordinates  is 

x  =  ( N  +  alt)  cos  <j)  cos  X  (1.21) 

y  =  (N  +  alt)  cos  (j)smX  (1-22) 

z  =  ((\-e2)N +  alt)sm(j)  (1-23) 

The  conversion  from  ECEF  to  ENU  coordinates  is  shown  below.  The 
subscript  obj  in  the  ECEF  coordinates  represents  the  object  of  interest  (in  our  case,  it  can 
be  the  satellite  or  GPS  receiver),  while  /  represents  the  local  reference  point. 

Est  =  -(xobj  -  x, )  sin  X  +  (yobj  -  yl )cosX  (1 .24) 

Nth  =  ~(xob/  -  x, )  sin  <f>cosX-  (yobj  -  y, )  sin  </>  sin  X  +  (zobj  -  z, )  cos  </>  (1.25) 

Up  =  (xobJ  - x, ) cos <f>cosX  +  (yobj  - y, ) cos </> sin T  +  (zobj  - z, ) sin </>  (1 .26) 
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B. 


DOP  CALCULATIONS 


1.  Identification  of  Visible  Satellites 

After  knowing  the  positions  of  all  the  satellites  at  the  given  time,  the  visible 
satellites  need  to  be  identified.  For  ease  of  calculation,  the  GPS  receiver's  position  is 
converted  from  the  geodetic  coordinates  (the  GPS  receiver's  position  is  entered  by  the 
user  in  geodetic  coordinates)  to  the  ECEF  coordinates  by  using  Equation  2.20  -  2.22. 

In  addition,  several  assumptions  are  made.  Firstly,  it  is  assumed  that  line  of  sight 
is  needed  for  the  satellite's  signal  to  be  visible  to  the  GPS  receiver.  Secondly,  the  Earth  is 
assumed  a  perfect  sphere.  This  differs  from  the  WGS  84  model  by  about  1  in  300  parts 
and  it  is  assessed  to  be  a  reasonable  assumption.  Thirdly,  the  GPS  receiver  is  assumed  to 
be  relatively  close  to  the  surface  of  the  Earth  such  that  the  field  of  view  of  the  sky  above 
the  receiver  is  constant  regardless  of  its  altitude. 

The  method  to  identify  visible  satellites  is  illustrated  in  Figure  3.  The  ECEF 
coordinates  of  the  GPS  receiver  and  satellite  are  represented  by  vectors  originating  from 

the  center  of  the  Earth,  OT  and  OS  respectively.  Any  obstruction  to  the  field  of  view  of 
the  entire  sky  above  the  receiver,  such  as  terrain,  is  represented  by  the  angle  /? .  a  ,  the 

angle  between  OT  and  TS  (where  TS  =  OS  -  OT  ),  is  given  by 

J  OT  TS  1 

a  =  cos  — - - r 

OT  TS 

V 

Hence,  the  satellite  is  visible  to  the  GPS  receiver  if  a  is  less  than  90°-  /3 . 


(1.27) 
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Satellite 


GPS  receiver 


a 


Obstruction 
to  Held  of  Mew 


Figure  3.  Schematic  for  Identifying  Visible  Satellite 


2.  DOP  Calculations 

With  the  visible  satellites  identified,  the  pseudoranges  between  the  GPS  receiver 
and  the  visible  satellites,  pi  are  calculated  as  shown  below.  The  subscript  i  represents  the 

numbering  of  each  visible  satellite,  while  r  represents  the  GPS  receiver.  Since  the  ENU 
frame  is  defined  from  the  local  reference  point,  Estr  and  Nth,-  are  zero,  while  Upr  is  the 
altitude  of  the  receiver. 


(1.28) 


The  directional  derivatives  of  East,  North,  Up  and  Time  in  the  ENU  frame  for 
each  satellite  are  then  calculated  as  shown  below. 


(1.29) 
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Ntht  -  Nthr 

^ Nth  ~ 

P, 

(1.30) 

UPi  -  Upr 

u»~  p, 

(1.31) 

A  ,=- 1 

(1.32) 

The  directional  derivatives  of  all  the  visible  satellites  fonn  a  (n  x  4)  matrix,  D, 
where  n  is  the  total  number  of  visible  satellites. 


^Estt 

D Nt\ 

DuPl 

D,  ' 

h 

D Est2 

DNth2 

DUp2 

Dtl 

(1.33) 

D Es,n 

D 'Nth„ 

DUP„ 

DJ 

Taking  the  inverse  of  DT  D  yields  a  (4  x  4)  matrix  where  the  diagonal  tenns 
gives  the  square  of  the  East  DOP  (XDOP),  North  DOP  (YDOP),  Vertical  DOP  (VDOP) 
and  Time  DOP  (TDOP)  as  shown  below  where  the  off-diagonal  tenns  are  not  shown  for 
clarity. 


(  XDOP 2 


V 


YDOP2 

VDOP2 

TDOP2  j 


(1.34) 


The  other  DOP  values,  namely  Horizontal  DOP  (HDOP),  Position  DOP  (PDOP) 
and  Geometric  DOP  (GDOP)  are  obtained  by  root  summing  the  appropriate  diagonal 
terms  in  ( DT  ■ D)  x  as  shown  below  [1]. 


HDOP  =  4  XDOP2  +  YDOP1  (1.35) 

PDOP  =  V  XDOP2  +  YDOP2  +  VDOP2  (1.36) 

GDOP  =  4  XDOP2  +  YDOP2  +  VDOP2  +  TDOP 2  (1.37) 
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III.  IMPLEMENTATION 


A.  APPROACH 

Matlab  was  selected  as  the  software  to  test  the  algorithm  of  the  DOP  calculator 
program  during  development  as  Matlab  is  relatively  easy  to  learn  and  manipulate 
compared  to  Visual  C++.  DOP  results  generated  by  the  algorithm  in  Matlab  were 
compared  with  DOP  calculations  generated  by  commercial  DOP  calculation  software  to 
check  the  algorithm. 

The  validated  algorithm  was  then  used  to  write  the  program  in  the  Visual  C++ 
environment  and  checked  again  with  the  results  generated  by  the  commercial  DOP 
calculation  software  to  ensure  that  the  algorithm  was  implemented  correctly  in  Visual 
C++.  In  addition,  the  option  to  plot  graphs  of  various  DOP  values  over  a  24-hour  period 
was  added  to  enhance  the  usability  of  the  program. 

B.  CODING  WITH  MATLAB  AND  VISUAL  C++ 

There  are  significant  differences  between  the  coding  language  in  Matlab  and 
Visual  C++.  For  example,  matrix  manipulation  is  integrated  into  the  Matlab  language, 
while  Visual  C++  is  more  generic  and  requires  additional  codes  to  be  written  in  order  to 
carry  out  matrix  operations.  The  same  is  true  for  graph  plotting. 

As  such,  instead  of  writing  C++  codes  for  matrix  manipulation  and  graph  plotting 
from  scratch,  open  source  C++  codes  available  from  the  Internet  are  modified  and 
incorporated  into  the  DOP  calculation  program.  The  Visual  C++  codes  for  DOP 
calculation  are  shown  in  Appendix  B  with  credits  given  to  the  parts  of  code  that  were 
modified  from  open  sources,  while  the  Matlab  codes  are  shown  in  Appendix  C. 

Figure  4  shows  the  screen  shot  of  the  Graphical  User  Interface  (GUI)  produced  by 
the  Visual  C++  codes.  The  user  inputs  are  grouped  into  "Time  Inputs"  and  "Position 
Inputs"  in  the  upper  part  of  the  window.  The  time  inputs  refer  to  the  local  time,  while  the 
position  inputs  refer  to  the  position  of  the  GPS  receiver  at  which  the  DOP  values  are  to 
be  calculated. 
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Pressing  the  "Update  Almanac"  button  will  download  the  latest  version  of  the 
almanac  file  from  http://www.navcen.uscg.gov/GPS/almanacs.htm,  which  is  the  United 
States  Coast  Guard  website,  to  the  folder  containing  the  rest  of  the  C++  codes.  The  codes 
stating  the  web  address  and  location  of  the  downloaded  file  are  in  the  C++  file 
"DOP  CalculatorDlg.cpp"  line  230  and  line  168  respectively. 

The  outputs  are  displayed  in  the  lower  part  of  the  window  after  the  "Compute 
DOP"  button  is  pressed.  The  left  portion  shows  the  various  calculated  DOP  values  for  the 
specified  time  inputs,  while  the  right  shows  the  graph  that  plots  the  DOP  values  over  a 
24-hour  period  on  the  specified  day.  By  default,  the  HDOP  is  plotted  but  other  DOP 
values  can  also  be  plotted  by  selecting  the  appropriate  radio  buttons. 


Time  Inputs 
Year (XXXX) 

Month  (1  to  12) 
Day  (1  to  31) 


2009 

2 

3 


Hour  (0  to  23)  0 _ 


Minute  (0  to  59)  0 


Timezone  (GMT/  Zulu  time  =  Ohr)  [  0  hr 
(EST  =  -5hr) 

Daylight  Saving  P 


Position  Inputs 

Longitude  (-180  to  180)  0 
(+East;  -West) 

deg 

Latitude  (-90  to  90)  0 

(+North;  -South) 

deg 

Altitude 


Obstruction  to 
Field  of  View 
(0  to  90) 


meters 

deg  from  the  horizontal 


Update  Almanac 


Outputs 


Compute  DOP 


XDOP  0.657 
YDOP  0.559 
VDOP  1.395 
TDOP  0.873 

HDOP  0.863 
PDOP  1.640 
GDOP  1.858 


Cancel 


Graph 

O  Plot  XDOP 
OPlot  YDOP 
©Plot  VDOP 
©  Plot  TDOP 

©Plot  HDOP 
O  Plot  PDOP 
©Plot  GDOP 


Variation  in  HDOP  over  24ht  Period 


Figure  4.  Graphical  User  Interface  of  the  DOP  Calculator 


27 


C.  CODE  TESTING 

DOP  values  generated  by  the  program  written  in  Matlab  and  Visual  C++  were 
compared  with  results  generated  by  the  commercial  software:  GPS  planning  software 
version  2.74  from  Trimble  Navigation  Limited,  which  is  publically  available  at  its 
website  [10].  Trimble’s  planning  software  was  chosen  for  a  couple  of  reasons. 

Firstly,  it  is  the  only  free  software  found  to  accept  the  same  input  requirements 
and  provide  direct  comparison  for  the  various  DOP  values  as  its  output.  Secondly, 
Trimble  is  a  listed  company  that  deals  with  GPS  applications,  giving  some  level  of 
credence  to  the  accuracy  of  the  results  generated  by  the  software. 

Only  GDOP,  PDOP,  HDOP,  TDOP  and  VDOP  were  compared,  as  these  are  the 
five  DOP  values  generated  by  the  Trimble  planning  software.  The  results  generated  by 
the  Trimble,  Matlab  and  Visual  C++  program  were  presented  on  graphs  and  tables  for 
comparison. 

Note  that  during  program  development,  results  generated  by  individual  sections  of 
the  codes,  such  as  the  calculation  of  the  position  of  satellite  in  the  ECEF  frame,  time 
conversion  from  Gregorian  to  Julian  date,  etc,  were  also  checked  against  external 
references  [1],  [6],  [11],  Hence,  the  combination  of  checking  individual  sections  and  the 
entire  program  with  different  external  references  ensures  that  the  DOP  calculation 
algorithm  was  implemented  correctly. 

1,  Matlab  Code  Testing 

Two  hundred  sixty  six  combinations  of  latitude  and  longitude  were  sampled  to 
test  the  codes  written  in  Matlab.  The  test  parameters  are  as  follows: 
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Table  3.  Test  Parameters  for  Matlab  Code 


Parameter 

Value 

Date  and  Time 

Feb  3,  2009,  OOOOhr  GMT 

Position  of  GPS  Receiver  in  Latitude 

-90°  (South)  to  90°  (North)  at  15°  interval 

Position  of  GPS  Receiver  in  Longitude 

-165°  (West)  to  180°  (East)  at  15°  interval 

Altitude  of  GPS  Receiver 

0  m 

Obstruction  to  Field  of  View,  /? 

jB  =  10° 

2.  Visual  C++ Code  Testing 

The  test  parameters  for  the  Visual  C++  codes  are  similar  to  those  for  the  Matlab 
codes,  except  that  there  were  less  data  points  as  extracting  DOP  values  from  Visual  C++ 
is  more  laborious.  62  combinations  of  latitude  and  longitude  were  sampled  and  the  test 
parameters  are  as  follows: 


Table  4.  Test  Parameters  for  Visual  C++  Code 


Parameter 

Value 

Date  and  Time 

Feb  3,  2009,  OOOOhr  GMT 

Position  of  GPS  Receiver  in  Latitude 

-90°  (South)  to  90°  (North)  at  30°  interval 

Position  of  GPS  Receiver  in  Longitude 

-150°  (West)  to  180°  (East)  at  30°  interval 

Altitude  of  GPS  Receiver 

0  m 

Obstruction  to  Field  of  View,  /? 

/?  =  10° 
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IV.  RESULTS 


A.  COMPARISON  OF  CALCULATED  DOP  VALUES 

1.  Matlab  vs  Trimble 

Of  the  266  combinations  of  latitude  and  longitude  sampled  around  the  Earth,  most 
of  the  differences  between  the  DOP  values  generated  by  the  Matlab  and  Trimble  program 
were  less  than  5%.  However,  39  of  the  samples  (about  15%  of  the  samples)  gave 
differences  that  were  between  5%  and  38%.  The  graphs  of  the  GDOP  values  are  shown  in 
Figure  5  below.  Although  there  were  differences  in  the  percentage  difference  for  various 
DOP  values,  GDOP  gives  a  good  representation  on  the  trend  since  it  is  derived  from  the 
rest  of  the  DOP  values.  The  complete  comparison  of  all  the  DOP  values  is  shown  in 
Appendix  D. 
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dOQD  dOQD  dOQD 


GDOP  vs  Longitude  (deg)  for  GDOP  vs  Longitude  (deg)  for 
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GDOP  vs  Longitude  (deg)  for 
Latitude  30  deg 


GDOP  vs  Longitude  (deg)  for 
Latitude  45  deg 
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Figure  5.  GDOP  generated  by  Matlab  and  Trimble  Program  at  Various  Positions  on 

Earth 


2.  Visual  C++  vs  Trimble 

From  the  62  samples,  most  of  the  differences  between  the  DOP  values  generated 
by  Visual  C++  and  Trimble  program  were  also  less  than  5%.  The  sample  points  that  gave 
more  than  5%  difference  were  the  same  as  those  between  the  Matlab  and  Trimble 
program.  This  is  expected  since  the  algorithm  in  the  Visual  C++  program  is  the  same  as 
the  one  in  Matlab  and  the  test  samples  for  the  Visual  C++  program  are  a  subset  of  those 
for  the  Matlab  program.  The  graphs  of  the  GDOP  values  are  shown  in  Figure  6  below, 
while  the  complete  comparison  of  all  the  DOP  values  is  shown  in  Appendix  E. 
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Figure  6.  GDOP  generated  by  Visual  C++  and  Trimble  Program  at  Various 

Positions  on  Earth 
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V.  ANALYSIS/  DISCUSSION  OF  RESULTS 


A.  COMPARISON  OF  CALCULATED  DOP  VALUES 

Since  the  results  generated  by  the  Visual  C++  program  are  a  subset  of  those 
generated  by  Matlab,  analysis  on  the  results  generated  by  the  Matlab  and  Trimble 
program  will  also  be  valid  for  the  Visual  C++  program.  Hence,  the  analysis  would  focus 
on  comparing  the  results  from  the  Matlab  and  Trimble  program  only. 

As  stated  in  the  previous  chapter,  for  the  DOP  values  generated  by  the  Matlab  and 
Trimble  program,  39  samples  had  differences  greater  than  5%.  It  was  found  that  the  GPS 
receiver  in  37  of  the  samples  had  different  number  of  satellites  in  its  field  of  view,  while 
two  samples  had  the  same  number  of  visible  satellites. 

1.  Difference  in  Number  of  Visible  Satellites 

In  the  37  samples,  all  of  them  had  the  number  of  visible  satellites  differed  by  one. 
In  every  case,  the  satellite  in  question  was  found  to  be  near  the  visibility  criterion.  The 
visibility  criterion  is  defined  by  the  obstruction  angle,  f3 .  In  this  case,  //  was  10°,  which 
meant  that  only  satellites  that  were  10°  above  the  horizon  were  considered  visible.  In 
other  words,  only  satellites  that  made  an  angle  of  less  than  80°  with  the  GPS  receiver 
(with  respect  to  the  Up-axis  in  ENU  frame;  see  angle  a  in  Figure  3)  were  in  the  field  of 
view  of  the  GPS  receiver.  Although  angle  a  ,  calculated  for  each  satellite  by  the  Matlab 
and  Trimble  program,  were  slightly  different  (differences  were  all  within  4°),  they  were 
enough  to  cause  one  satellite  to  be  included  or  excluded  from  the  field  of  view  of  the 
GPS  receiver  in  these  samples. 

a.  Leap  Seconds 

It  was  found  that  in  one  of  the  37  samples,  this  discrepancy  could  be 
accounted  for  by  taking  the  effects  of  leap  seconds  into  account.  When  two  leap  seconds 
were  removed  from  the  Matlab  calculations  (in  the  time  conversion  algorithm  to  convert 
to  GPS  time),  the  difference  in  the  number  of  visible  satellites  of  that  sample  (Latitude  - 
90°  (South),  Longitude  0°)  disappeared.  For  this  sample  point,  satellite  number  1 1  went 
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from  visible  to  invisible  in  the  Matlab  calculation  and  this  changed  the  DOP  values  such 
that  the  differences  compared  to  the  results  from  the  Trimble  program  decreased  from 
about  24%  to  less  than  4%.  As  the  source  code  for  the  Trimble  program  is  not  available, 
it  is  speculated  that  this  version  of  the  Trimble  program  was  developed  before  2005  and 
as  a  result,  the  two  additional  leap  seconds,  which  were  added  on  Jan  1,  2005  and  Jan  1, 
2009  [12],  were  not  accounted  for. 

b.  Visibility  Algorithm 

For  the  other  36  samples,  when  the  number  of  visible  satellites  were 
deliberately  made  to  be  the  same  as  those  seen  by  the  Trimble  program  by  adjusting  the 
obstruction  angle,  ft ,  the  differences  were  found  to  drop  to  less  than  5%.  Hence,  the 
discrepancies  in  these  samples  could  be  attributed  to  the  difference  in  the  algorithm  used 
by  the  Trimble  and  Matlab  program  to  detennine  the  visibility  of  satellites.  It  is  possible 
it  is  due  to  different  assumptions  made  on  the  visibility  algorithm. 

Interestingly,  it  was  found  that  one  of  the  samples  (Latitude  60°  (North), 
Longitude  90°  (East):  not  among  the  37  samples  mentioned  earlier)  had  a  difference  of 
less  2%  in  the  DOP  values  although  the  number  of  visible  satellites  also  differed  by  one. 
On  closer  examination,  it  was  found  that  there  were  two  satellites  (satellite  number  15 
and  26)  close  to  each  other  that  were  near  the  visibility  criterion.  Since  DOP  is  dependent 
on  the  angular  separation  of  the  visible  satellites,  the  loss  of  one  of  these  two  satellites 
from  the  field  of  view  caused  little  change  to  the  DOP.  In  this  case,  the  angular  separation 
between  the  two  satellites  was  about  18°.  Satellite  number  26  was  in  the  field  of  view  in 
the  Trimble  program,  while  it  was  just  out  of  it  in  the  Matlab  program. 

Details  on  the  comparison  of  the  angles  of  the  visible  satellites  are  shown 
in  Appendix  F.  Note  that  the  angles  that  caused  the  discrepancies  are  shaded  in  grey  for 
clarity. 
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2.  No  Difference  in  Number  of  Visible  Satellites 

Having  accounted  for  the  large  DOP  percentage  differences  in  37  of  the  39 
samples,  the  remaining  two  samples  (Latitude  -45°  (South),  Longitude  60°  (East)  & 
Latitude  0°,  Longitude  90°  (East))  had  the  same  number  of  visible  satellites  from  both 
programs. 

On  closer  examination,  it  was  found  that  although  the  total  number  of  visible 
satellites  was  the  same,  the  actual  satellites  seen  by  the  GPS  receiver  were  different.  By 
sheer  coincidence,  one  of  the  satellites  was  considered  visible  in  the  Trimble  program  but 
not  in  Matlab  program,  while  another  satellite  was  considered  visible  in  the  Matlab 
program  but  not  in  the  Trimble  program. 

For  the  sample  point  Latitude  -45°  (South),  Longitude  60°  (East),  satellite  number 
26  was  considered  visible  by  the  Trimble  program  but  was  marginally  out  of  the  visibility 
criterion  in  the  Matlab  program,  while  the  reverse  is  true  for  satellite  number  23.  This 
caused  the  DOP  values  to  differ  by  about  15%. 

For  Latitude  0°,  Longitude  90°  (East),  satellite  number  1 1  was  considered  visible 
by  the  Trimble  program  but  invisible  by  the  Matlab  program,  while  the  reverse  is  true  for 
satellite  number  26.  This  resulted  in  about  7%  difference  in  DOP  values. 

3.  Effect  of  One  Visible  Satellite  on  DOP 

The  effect  of  an  addition  or  loss  of  a  visible  satellite  (or  different  visible  satellite 
for  the  case  where  the  total  number  of  visible  satellites  is  the  same)  on  the  DOP  value  is 
dependent  on  the  angular  separation  between  that  satellite  and  the  rest  of  the  visible 
satellites.  If  none  of  the  remaining  visible  satellites  were  near  that  satellite,  the  effect  of 
seeing  or  losing  sight  of  the  additional  satellite  near  the  visibility  criterion  would  make  a 
disproportionally  big  difference  to  the  DOP  values.  The  greatest  difference  in  the  DOP 
values  was  found  to  be  about  38%  due  to  a  loss  of  one  visible  satellite  (Latitude  -60° 
(South),  Longitude  75°  (North)  &  Latitude  -60°  (South),  105°  (East),  where  there  were  9 
other  visible  satellites  remaining  in  the  field  of  view). 
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Conversely,  if  there  were  visible  satellites  near  that  particular  satellite,  the  effect 
of  seeing  or  losing  sight  of  the  additional  satellite  would  be  small.  This  is  evident  in  the 
sample  point  Latitude  60°  (North),  Longitude  90°  (East)  having  a  difference  of  less  than 
2%  in  the  DOP  values  although  the  number  of  visible  satellite  also  differed  by  one. 

B.  EFFECT  OF  USING  OUTDATED  ALMANAC  DATA 

The  almanac  file  is  updated  almost  daily.  However,  there  might  be  circumstances 
where  the  most  updated  almanac  file  could  not  be  accessed.  As  such,  the  effect  of  using 
outdated  almanac  data  on  DOP  values  was  studied. 

1.  Approach 

In  the  study,  an  arbitrary  date  was  chosen  and  the  almanac  file  was  used  as  the 
benchmark.  The  chosen  date  in  this  study  was  Jul  28,  2008  (Almanac  file:  Week  466, 
Time  of  Applicability  319488).  Using  this  almanac  file,  the  DOP  values  on  dates  that 
were  1,7,  14  and  30  days  later  were  calculated.  The  position  of  the  GPS  receiver  was 
fixed  at  Latitude  0°,  Longitude  0°  at  altitude  0  m  and  the  obstruction  angle,  ft  was  set  at 
10°.  These  DOP  values  were  then  compared  with  the  values  obtained  from  using  the 
correct  almanac  files  for  the  respective  dates. 

2.  No  Difference  in  Number  of  Visible  Satellites 

It  was  found  that  using  outdated  almanac  data  had  little  effect  on  the  DOP  values 
if  the  total  number  of  visible  satellites  remained  the  same.  Even  when  the  almanac  was 
30-days  old,  the  maximum  difference  for  these  sample  points  was  less  than  3%. 
However,  there  was  a  general  trend  that  the  differences  became  larger  as  the  almanac 
data  became  more  outdated.  This  is  expected  as  the  more  outdated  almanac  data  would  be 
expected  to  have  less  accurate  corrections  on  the  orbital  deviations  of  the  satellites 
compared  to  more  recent  almanac  data. 

3.  Difference  in  Number  of  Visible  Satellites 

However,  in  the  case  where  the  total  number  of  visible  satellites  is  different,  large 
errors  may  be  possible.  In  this  case,  the  differences  in  DOP  values  were  found  to  be  as 


39 


high  as  25%  in  one  of  the  sample  points  when  the  30-days  old  almanac  was  used. 
Previous  sections  in  the  report  had  already  discussed  the  effect  of  having  one  more  or  less 
visible  satellite  on  the  DOP. 

However,  note  that  the  reason  for  the  difference  in  the  total  number  of  visible 
satellite  is  different  when  it  comes  to  the  use  of  outdated  almanac  data.  In  this  case,  the 
use  of  such  data  in  the  calculations  resulted  in  satellite  positions  that  were  slightly 
different  compared  to  the  satellite  positions  if  the  correct  almanac  data  were  used.  As  a 
result,  when  a  satellite  is  near  the  visibility  criterion,  the  exact  time  when  the  satellite 
goes  into  or  out  of  the  field  of  view  is  different  depending  on  the  almanac  data  used.  The 
difference  in  the  exact  time  increases  when  the  almanac  is  more  outdated  hence,  the 
probability  of  getting  the  wrong  set  of  satellites  in  the  field  of  view  increases  as  the 
almanac  data  becomes  more  outdated.  Hence,  wherever  possible,  the  most  recent  almanac 
data  should  be  used. 

Details  on  the  DOP  values  calculated  from  outdated  almanac  data  are  shown  in 
Appendix  G. 

C.  AVERAGE  HDOP  AND  VDOP 

In  the  event  that  time  and  location  cannot  be  specified  ahead  of  time,  it  is  still 
useful  to  have  a  "rule  of  thumb"  on  the  DOP  values  during  mission  planning.  The  most 
useful  DOP  values  for  GPS-guided  weapons  are  the  HDOP  and  VDOP  as  these  values 
are  used  in  calculations  to  estimate  the  Circular  and  Height  Error  Probable  (CEP  &  HEP). 
The  current  "rule  of  thumb"  regarding  the  ratio  between  VDOP  and  HDOP  is  that  the 
VDOP  is  twice  the  value  of  HDOP  [13]. 

1.  Approach 

In  order  to  verify  the  current  "rule  of  thumb"  as  well  as  to  establish  another  for  the 
global  average  value  of  HDOP  and  VDOP,  one  million  iterations  of  the  DOP  calculation 
were  done  to  obtain  the  average  DOP  values  for  random  positions  on  the  Earth's  surface 
(i.e.,  altitude  is  zero)  at  random  time  for  a  specified  day.  Six  different  days  were  sampled 
and  the  mean  of  the  average  DOP  values  from  each  day  is  taken  to  be  the  global  average. 
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Three  values  of  obstruction  angle,  (3  =  0°,  10°  and  15°,  were  used  to  get  three  sets 
of  values  for  HDOP  and  VDOP.  /?  =  0°  represented  the  most  optimistic  case  where  there 
is  no  obstruction  to  the  field  of  view  of  the  GPS  receiver.  This,  however,  may  not  be 
realistic  in  many  situations  hence,  f3  =  10°  and  15°  were  used  to  represent  cases  that  are 
more  realistic. 

In  order  to  get  the  statistical  spread  of  the  ratio  between  VDOP  and  HDOP,  the 
VDOP  and  HDOP  generated  by  the  Matlab  program  that  were  used  to  compare  against 
results  generated  by  the  Trimble  program  in  the  previous  sections  were  extracted  and  the 
ratio  between  the  two  DOP  values  were  calculated  for  each  sample  point. 

2.  Results 

The  results  on  the  global  average  VDOP  and  HDOP  are  shown  in  Table  5  below. 


Table  5.  Global  Average  VDOP  and  HDOP 


p 

Global  Average 

VDOP 

Global  Average 

HDOP 

Ratio  between  Global  Average 

VDOP  and  Global  Average  HDOP 

0° 

1.2 

0.8 

1.5  :  1 

10° 

1.8 

1.0 

1.8  :  1 

15° 

2.4 

1.2 

2.0  :  1 

From  the  results  above,  the  "rule  of  thumb"  that  the  VDOP  value  is  twice  the 
value  of  HDOP  was  found  to  be  most  accurate  when  (3  is  15°.  The  ratio  between  VDOP 
and  HDOP  was  observed  to  increase  as  f3  increase.  However,  note  that  the  trend  of 
increasing  ratio  would  not  continue  for  large  [3 .  This  is  because  fewer  satellites  would  be 
visible  as  f3  increases  and  there  would  come  a  point  where  there  were  less  than  four 
satellites  in  the  field  of  view-  the  minimum  number  to  calculate  the  GPS  receiver's 
position.  Details  on  the  global  average  of  DOP  are  shown  in  Appendix  H. 
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When  the  extracted  VDOP  and  HDOP  values  from  the  266  samples  generated  by 
the  Matlab  program  were  examined,  the  ratio  between  VDOP  and  HDOP  for  /3  =  10° 
was  found  to  vary  widely  from  about  0.8  to  3.6  with  an  average  and  standard  deviation  of 
about  1.8  and  0.4  respectively. 

Although  the  "rule  of  thumb"  that  the  VDOP  is  twice  the  value  of  HDOP 
remained  to  be  a  simple  and  reasonable  guide,  caution  should  be  exercised  as  the  spread 
of  the  ratio  was  shown  to  be  fairly  large.  The  study  also  gives  planners  the  flexibility  to 
use  other  ratios  if  more  information  is  known  about  the  obstructions  to  the  field  of  view 
of  the  GPS  receiver  in  the  operating  area.  Moreover,  the  study  provided  planners  with  the 
average  values  of  VDOP  and  HDOP  that  can  be  used  if  actual  DOP  values  are  not  readily 
available. 

Details  on  the  results  are  shown  in  Appendix  I.  Note  that  the  highest  and  lowest 
ratio  between  VDOP  and  HDOP  are  shaded  in  grey  for  clarity. 


42 


THIS  PAGE  INTENTIONALLY  LEFT  BLANK 


43 


VI.  CONCLUSION 


A  DOP  calculation  program  that  extracts  satellite  infonnation  from  the  SEM 
almanac  file  had  been  successfully  developed  in  the  Visual  C++  environment. 
Comparison  of  the  results  with  a  commercial  DOP  calculation  program,  the  Trimble 
Planning  software,  showed  that  the  algorithm  had  been  implemented  correctly. 

The  effect  of  a  difference  in  the  number  of  visible  satellites  on  DOP  had  been 
investigated.  It  was  found  the  addition  or  loss  of  one  satellite  from  the  field  of  view  has 
the  potential  to  result  in  disproportional  effects  on  the  DOP  values  depending  on  the 
configuration  of  the  rest  of  the  visible  satellites.  The  change  in  DOP  values  had  been 
found  to  be  as  large  as  38%. 

The  effects  of  using  outdated  almanac  data,  up  to  30-days  old,  had  also  been 
studied.  It  was  found  that  the  effect  on  DOP  is  small  unless  the  number  of  satellites  in  the 
field  of  view  of  the  GPS  receiver  is  different.  The  more  outdated  the  almanac  data,  the 
higher  the  probability  of  encountering  cases  of  different  number  of  visible  satellites. 

It  should  be  highlighted  that,  although  discrepancies  in  number  of  visible  satellites 
could  be  attributed  to  small  differences  in  the  algorithm  between  the  Matlab  and  Trimble 
program,  the  level  of  uncertainty  during  actual  GPS  usage  can  be  expected  to  have  a 
much  larger  effect.  Isolated  buildings  or  hills  may  block  the  signal  from  an  otherwise 
visible  satellite,  causing  the  effective  number  of  visible  satellites  to  decrease.  Hence,  an 
important  take-away  is  to  realize  the  potential  uncertainty  in  the  DOP  values  where 
conditions  on  the  ground  can  deviate  from  the  theoretical  ideal  calculated  in  the  program. 

Finally,  the  relationship  between  HDOP  and  VDOP  had  been  studied.  A  "rule  of 
thumb"  on  the  global  average  values  of  HDOP  and  VDOP  was  established.  This  would 
allow  planners  to  have  a  better  estimate  on  weapon  delivery  accuracy  when  specific  time 
and  position  of  target  are  unknown. 
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A.  FUTURE  WORK 

1.  Integration  of  Code  into  JWS 

The  next  stage  of  development  of  the  program  is  to  integrate  the  codes  for  DOP 
calculation  into  the  JWS.  As  the  JWS  operates  in  a  secure  network,  direct  access  to 
Internet  is  not  possible.  As  such,  the  web  address  to  download  the  latest  almanac  file  as 
well  as  the  location  to  retrieve  the  stored  almanac  file  in  the  secure  network  would  need 
to  be  updated  in  the  code. 

2.  Enhancement  of  GUI 

Throughout  the  program  development,  the  user-friendliness  of  the  GUI  had  been  a 
consideration.  However,  due  to  limited  time,  feedback  from  potential  users  of  the 
program  had  not  been  solicited.  Future  versions  of  the  program  should  elicit  feedback 
from  users  to  improve  the  usability  and  functionality  of  the  program. 
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APPENDIX  A.  EXAMPLES  OF  ALMANAC  FILE 


A.  EXAMPLE  OF  ALMANAC  FILE  IN  SEM  (.AL3)  FORMAT 

31  CURRENT.  ALM 
493  405504 

2 

61 

0 

8. 9578628540039 IE-03  -3.29971313476562E-04  -2.5429471861571  IE-09 
5. 1 5354736328 125E+03  5.08018851280212E-01  8.63925099372864E-01 
-3.86929750442505E-01  1.58309936523438E-04  0.00000000000000E+00 
0 
9 

3 

33 
0 

1 .16286277770996E-02  -5.1 0406494 140625E-03  -2.582964953035 12E-09 
5. 1 53666992 18750E+03  1.43044233322144E-01  2.83277988433838E-01 
-5.01049160957336E-01  3.49044799804688E-04  3.63797880709171E-12 
0 
9 

4 

34 
0 

8.6336 1358642578E-03  -5.72204589843750E-04  -2.561 13708019257E-09 
5.15361083984375E+03  5.13915061950684E-01  1.50653243064880E-01 
5. 1 579928398 1323E-01  -3.26156616210938E-04  - 1 .455 1 9 1 52283669E- 1 1 
0 
9 
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5 


35 

1 

9.39798355 102539E-03  8.81 195068359375E-04  -2.50292941927910E-09 
5. 1 53577 1 4843 750E+03  -1.85006141662598E-01  4. 12206888 1 98853E-0 1 
3.9079093933 1055E-01  -6.39915466308594E-04  -4.3655745685 1 006E- 1 1 
0 
9 


6 

36 

0 

5.70201873779297E-03 
5.15364599609375E+03 
3. 0269 1 57886505  IE-01 
0 
9 


-2.725601 19628906E-03  -2.561 13708019257E-09 
1. 6426539421 08 15E-01  -4.65690255165100E-01 
6. 77 1 0876464843 8E-05  - 1 .455 1 9152283 669E- 1 1 


7 

48 

0 

2.29692459106445E-03 

5.15363134765625E+03 

-9.65989708900452E-01 

0 

10 


7.37380981445312E-03  -2.57568899542093E-09 
-4.88172769546509E-01  9.30595278739929E-01 
2.28881835937500E-05  0.00000000000000E+00 


8 

38 

0 

1 .047 1 82083 12988E-02  1.41 143798828125E-02  -2.49929144047201E-09 
5.15364062500000E+03  -4.69682693481445E-01  9.47455763816833E-01 
8.207 14473 724365E-01  -1.89781 188964844E-04  0.00000000000000E+00 
0 
9 
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9 


39 

0 

2.0 175933  83  78906E-02 
5. 1 53604492 18750E+03 
7. 3905277252 1973E-01 
0 
9 


1 .02252960205078E-02  -2.550223 14377 129E-09 
-4.99263 1673  8 1287E-01  4.6454501 1520386E-01 
3.91006469726562E-05  3 .63 797880709 1 7 IE- 12 


10 

40 

0 

8.52298736572266E-03  5.58662414550781E-03  -2.46654963 120818E-09 
5. 1 5373 876953 125E+03  8.5304892063 1409E-01  1.88068747520447E-01 
-1 . 12927794456482E-01  -1.04904 1 74804688E-05  0.00000000000000E+00 
0 
9 


11 

46 

0 

9.18340682983398E-03 
5. 1 53588867 18750E+03 
9. 90043640 1 367 19E-01 
0 
9 


-1.68552398681 64 IE-02  -2.74667399935424E-09 
4.460672 140 12 146E-01  2.17745780944824E-01 
8.58306884765625E-06  0.00000000000000E+00 


12 

58 

0 

3.19671630859375E-03 
5.15368066406250E+03 
-9.0521 33560 18066E-01 
0 
10 


7.56835937500000E-03  -2.42289388552308E-09 
-1.5343201 1604309E-01  -2.5871 1695671082E-01 
-3. 18527221 679688E-04  3.63797880709171E-12 
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13 


43 

0 

3. 744 1253662 1094E-03 
5. 1 537 1 533203 125E+03 
-1.703 164577484 13E-01 
0 
9 


1 .68495 1 78222656E-02  -2.4920 1 548285782E-09 
-8.04313659667969E-01  4.79433 17890 1672E-01 
2.88009643554688E-04  0.00000000000000E+00 


14 

41 

0 

4.41741943359375E-03 

5.15352929687500E+03 

-2.97486662864685E-01 

0 

9 


1 .4808654785 1562E-02  -2.51020537689328E-09 
-8.1030011 1770630E-01  -6.68291091918945E-01 
-1. 5830993652343 8E-04  3.63797880709171E-12 


15 

55 

0 

1.40905380249023E-03 
5. 1 535283203 1250E+03 
-1.366590261 45935E-01 
0 
10 


4. 673004 15039062E-03  -2.601 15484707057E-09 
-8.25886249542236E-01  -2.36185193061829E-01 
-2. 52723693 847656E-04  -3.63797880709171E-12 


16 

56 

0 

5.04970550537 109E-03 
5. 1 537 1 533203 125E+03 
2. 9290831 0890 198E-01 
0 
9 


8.05664062500000E-03  -2.42653186433017E-09 
-1. 4809322357 1777E-01  -1.22860670089722E-01 
7.72476 196289062E-05  -3.63797880709171E-12 
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17 


53 

0 

4. 197 1206665039 IE-03 
5. 1 535786 1 328 125E+03 
-5. 2773  83327484 13E-02 
0 
10 


5.53703308 105469E-03  -2.48110154643655E-09 
1.83424353599548E-01  -8.75360608100891E-01 
4. 5776367 1875000E-05  0.00000000000000E+00 


18 

54 

0 

9.91773605346680E-03  5.1 1 169433593750E-04  -2.51020537689328E-09 
5.15370654296875E+03  8.553 1 6 1 62 1093 75E-0 1  -7.76121497154236E-01 
3.33639025688 17 IE-01  -8.58306884765625E-05  3.63797880709171E-12 
0 
9 


19 

59 

0 

5.12599945068359E-03 

5.15372167968750E+03 

-2.88856863975525E-01 

0 

9 


4.80270385742188E-03  -2.49565346166492E-09 
2.0136451721 1914E-01  -1. 245497465 13367E-01 
3. 24249267578 125E-05  0.00000000000000E+00 


20 

51 

0 

3.94058227539062E-03 
5.15363427734375E+03 
3. 8049054 1 458 130E-01 
0 
9 


3. 7 1 93298339843 8E-04  -2.5 1748 133450747E-09 
8.38206768035889E-01  4. 18452620506287E-01 
9.1 5527343 750000E-05  0.00000000000000E+00 


50 


21 

45 

0 

1 .46293640 1 367 19E-02  -2.63404846191406E-03 
5.15357958984375E+03  5. 1 833677291 870 1E-0 1 
7. 3988 16347 1221 9E-01  3.05175781250000E-05 
0 
9 


22 

47 

1 

4. 88758087 158203E-03 
5. 1 536 1 767578 125E+03 
-2.4381 1607360840E-02 
0 
9 

23 
60 
0 

5.77783584594727E-03 
5.15360058593750E+03 
-4. 604705572 12830E-01 
0 
9 

24 
24 
0 

6.956 100463  867 19E-03 
5.15360839843750E+03 
3.0127465724945  IE-01 
0 
9 


-2.32696533203 125E-04 
8.57203364372253E-01 
2.05993652343750E-04 


-2.57568899542093E-09 
-8.43 170762062073E-01 
0.00000000000000E+00 


-2.5 13  8433557003  7E-09 
-5.80405235290527E-01 
0.00000000000000E+00 


9.09423828125000E-03  -2.56477505899966E-09 
-8.2002735 1379395E-01  9.17497277259827E-01 
3. 8623 8098 14453 IE-04  0.00000000000000E+00 


2.66265869 140625E-03  -2.52475729212165E-09 
5. 2797 1267700 195E-01  -2.03566551208496E-01 
1. 63078308 105469E-04  3.63797880709171E-12 
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25 

25 

1 


1.1797428131 1035E-02  8.1 19583 12988281E-03  -2.57568899542093E-09 
5. 1 5362548828 125E+03  -5.18336057662964E-01  -3.91333818435669E-01 
4.67219591 140747E-01  2.23159790039062E-04  2.54658516496420E-11 
0 
9 


26 

26 

0 

1.95465087890625E-02 
5. 1 5360 10742 1875E+03 
-6.00259423255920E-01 
0 
9 


1 .6 1 743 1 64062500E-02  -2.4920 1 548285782E-09 
-8.04983854293823E-01  3.08019638061523E-01 
4. 196 1669921 8750E-05  3 .63 797880709 1 7 IE-12 


27 

27 

1 

2.10742950439453E-02 
5.15513378906250E+03 
1.71 7686653 1 3  72 1E-0 1 
0 
9 


9.5996856689453  IE-03  -2.53567 122854292E-09 
-5.06914854049683E-01  -5.25584697723389E-01 
2.1 934509277343 8E-05  3.63797880709171E-12 


28 

44 

1 


1 .4636993408203  IE-02  7.48062133789062E-03  -2.42289388552308E-09 
5.15359619140625E+03  -1.45439743995667E-01  -6.54806971549988E-01 
8.22650 194 1 6809 IE-02  -2.38418579101562E-05  0.00000000000000E+00 
0 
9 


52 


29 


57 

0 

3.26204299926758E-03  5.59425354003906E-03  -2.48837750405073E-09 
5. 1 5368408203 125E+03  1 .8602883  8 1 57654E-0 1  -4.27628 1 595230 1 OE-0 1 
7.56655812263489E-01  -1.23977661 132812E-05  3.63797880709171E-12 
0 
10 


30 

30 

0 

1.1 1746788024902E-02 
5.15359765625000E+03 
2. 1 978664398 1934E-01 
0 
9 


2. 7790069580078 IE-03  -2.48110154643655E-09 
- 1.6805875301 36 11E-01  4.49827551841736E-01 
1.24931 33544921 9E-04  3.63797880709171E-12 


31 

52 

1 

7. 1 887969970703  IE-03 
5.15359130859375E+03 
-9.893 19086074829E-01 
0 
10 


9. 09423828 125000E-03  -2.553861 12257838E-09 
-4.89016652107239E-01  -3.80007505416870E-01 
-5. 054473 876953 12E-05  0.00000000000000E+00 


32 

23 

0 

1.3553 1425476074E-02 
5.15360058593750E+03 
-7. 171 325683593  75E-01 
0 
9 


7.39669799804688E-03  -2.44108377955854E-09 
8.723 16479682922E-01  -3.82601380348206E-01 
2. 975463 8671 8750E-04  -3.63797880709171E-12 
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B.  EXAMPLE  OF  ALMANAC  FILE  IN  YUMA  (.ALM)  FORMAT 

********  Week  493  almanac  for  PRN-02  ******** 

ID:  02 

Health:  000 

Eccentricity:  0.8957862854E-002 

Time  of  Applicability(s):  405504.0000 
Orbital  Inclination(rad):  0.9414405823 
Rate  of  Right  Ascen(r/s):  -0.7992639439E-008 
SQRT(A)  (m  1/2):  5153.547363 

Right  Ascen  at  Week(rad):  0.1595988274E+001 
Argument  of  Perigee(rad):  2.714100718 
Mean  Anom(rad):  -0.1215575695E+001 

Af0(s):  0.1583099365E-003 

Afl(s/s):  0.0000000000E+000 

week:  493 

********  Week  493  almanac  for  PRN-03  ******** 

ID:  03 

Health:  000 

Eccentricity:  0.1 162862778E-001 

Time  of  Applicability(s):  405504.0000 
Orbital  Inclination(rad):  0.9264431000 
Rate  of  Right  Ascen(r/s):  -0.81 16330719E-008 
SQRT(A)  (m  1/2):  5153.666992 

Right  Ascen  at  Week(rad):  0.4493867159E+000 
Argument  of  Perigee(rad):  0.889944077 
Mean  Anom(rad):  -0.1574092388E+001 

Af0(s):  0.3490447998E-003 

Afl(s/s):  0.3637978807E-01 1 

week:  493 

********  Week  493  almanac  for  PRN-04  ******** 

ID:  04 

Health:  000 

Eccentricity:  0.8633613586E-002 

Time  of  Applicability(s):  405504.0000 
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Orbital  Inclination(rad):  0.9406795502 
Rate  of  Right  Ascen(r/s):  -0.8043571 142E-008 
SQRT(A)  (m  1/2):  5153.610840 

Right  Ascen  at  Week(rad):  0.1614511 728E+00 1 
Argument  of  Perigee(rad):  0.473291 159 
Mean  Anom(rad):  0.1620431185E+001 

Af0(s):  -0.3261 566 162E-003 

Afl(s/s):  -0.1455191523E-010 

week:  493 

********  Week  493  almanac  for  PRN-05  ******** 

ID:  05 

Health:  000 

Eccentricity:  0.939798355  IE-002 

Time  of  Applicability(s):  405504.0000 

Orbital  Inclination(rad):  0.9452457428 

Rate  of  Right  Ascen(r/s):  -0.7858034223E-008 

SQRT(A)  (m  1/2):  5153.577148 

Right  Ascen  at  Week(rad):  -0.581213951 1E+000 

Argument  of  Perigee(rad):  1.294986129 

Mean  Anom(rad):  0.1227705956E+001 

Af0(s):  -0.6399154663E-003 

Afl(s/s):  -0.4365574569E-010 

week:  493 

********  Week  493  almanac  for  PRN-06  ******** 

ID:  06 

Health:  000 

Eccentricity:  0.5702018738E-002 

Time  of  Applicability(s):  405504.0000 

Orbital  Inclination(rad):  0.9339141846 

Rate  of  Right  Ascen(r/s):  -0.8050847100E-008 

SQRT(A)  (m  1/2):  5153.645996 

Right  Ascen  at  Week(rad):  0.5160549879E+000 

Argument  of  Perigee(rad):  -1.463009119 

Mean  Anom(rad):  0.9509336948E+000 
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AfO(s): 

Afl(s/s): 

week: 


0.6771 087646E-004 


-0.1455191523E-010 
493 

********  Week  493  almanac  for  PRN-07  ******** 

ID:  07 

Health:  000 

Eccentricity:  0.2296924591E-002 

Time  of  Applicability(s):  405504.0000 

Orbital  Inclination(rad):  0.9656429291 

Rate  of  Right  Ascen(r/s):  -0.8090864867E-008 

SQRT(A)  (m  1/2):  5153.631348 

Right  Ascen  at  Week(rad):  -0.1533640027E+001 

Argument  of  Perigee(rad):  2.923551321 

Mean  Anom(rad):  -0.3034746170E+001 

Af0(s):  0.228881 8359E-004 

Afl(s/s):  0.0000000000E+000 

week:  493 

********  Week  493  almanac  for  PRN-08  ******** 

ID:  08 

Health:  000 

Eccentricity:  0.1047182083E-001 

Time  of  Applicability(s):  405504.0000 

Orbital  Inclination(rad):  0.9868202209 

Rate  of  Right  Ascen(r/s):  -0.7847120287E-008 

SQRT(A)  (m  1/2):  5153.640625 

Right  Ascen  at  Week(rad):  -0.1475551724E+001 

Argument  of  Perigee(rad):  2.976520061 

Mean  Anom(rad):  0.2578350544E+001 

Af0(s):  -0.189781 1890E-003 

Afl(s/s):  0.0000000000E+000 

week:  493 

********  Week  493  almanac  for  PRN-09  ******** 
ID:  09 
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Health:  000 

Eccentricity:  0.2017593384E-001 

Time  of  Applicability(s):  405504.0000 
Orbital  Inclination(rad):  0.9746017456 
Rate  of  Right  Ascen(r/s):  -0.80144673 12E-008 
SQRT(A)  (m  1/2):  5153.604492 

Right  Ascen  at  Week(rad):  -0.1568481445E+001 
Argument  of  Perigee(rad):  1.45941 1 144 
Mean  Anom(rad):  0.2321802735E+001 

Af0(s):  0.3910064697E-004 

Afl(s/s):  0.3637978807E-01 1 

week:  493 

********  Week  493  almanac  for  PRN-10  ******** 

ID:  10 

Health:  000 

Eccentricity:  0.8522987366E-002 

Time  of  Applicability(s):  405504.0000 

Orbital  Inclination(rad):  0.9600296021 

Rate  of  Right  Ascen(r/s):  -0.7745256880E-008 

SQRT(A)  (m  1/2):  5153.738770 

Right  Ascen  at  Week(rad):  0.2679932237E+001 

Argument  of  Perigee(rad):  0.590835452 

Mean  Anom(rad):  -0.354773 1638E+000 

Af0(s):  -0.1 04904 1748E-004 

Afl(s/s):  0.0000000000E+000 

week:  493 

********  Week  493  almanac  for  PRN- 11  ******** 

ID:  11 

Health:  000 

Eccentricity:  0.9183406830E-002 

Time  of  Applicability(s):  405504.0000 

Orbital  Inclination(rad):  0.8895263672 

Rate  of  Right  Ascen(r/s):  -0.8625647752E-008 

SQRT(A)  (m  1/2):  5153.588867 
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Right  Ascen  at  Week(rad):  0.1401361465E+001 
Argument  of  Perigee(rad):  0.684068561 
Mean  Anom(rad):  0.3 1 1 03 1 3773E+00 1 

AfO(s):  0.8583068848E-005 

Afl(s/s):  0.0000000000E+000 

week:  493 

********  Week  493  almanac  for  PRN-12  ******** 

ID:  12 

Health:  000 

Eccentricity:  0.31 967 16309E-002 

Time  of  Applicability(s):  405504.0000 

Orbital  Inclination(rad):  0.9662551880 

Rate  of  Right  Ascen(r/s):  -0.7617927622E-008 

SQRT(A)  (m  1/2):  5153.680664 

Right  Ascen  at  Week(rad):  -0.4820208549E+000 

Argument  of  Perigee(rad):  -0.812766790 

Mean  Anom(rad):  -0.284381 163 1E+001 

Af0(s):  -0.3 1852722 17E-003 

Afl(s/s):  0.3637978807E-011 

week:  493 

********  Week  493  almanac  for  PRN-13  ******** 

ID:  13 

Health:  000 

Eccentricity:  0.3744125366E-002 

Time  of  Applicability(s):  405504.0000 

Orbital  Inclination(rad):  0.9954128265 

Rate  of  Right  Ascen(r/s):  -0.7832568372E-008 

SQRT(A)  (m  1/2):  5153.715332 

Right  Ascen  at  Week(rad):  -0.2526825905E+001 

Argument  of  Perigee(rad):  1.506183743 

Mean  Anom(rad):  -0.5350649357E+000 

Af0(s):  0.2880096436E-003 

Afl(s/s):  0.0000000000E+000 

week:  493 
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********  Week  493  almanac  for  PRN-14  ******** 

ID:  14 

Health:  000 

Eccentricity:  0.441 74 19434E-002 

Time  of  Applicability(s):  405504.0000 

Orbital  Inclination(rad):  0.9890003204 

Rate  of  Right  Ascen(r/s):  -0.7883500075E-008 

SQRT(A)  (m  1/2):  5153.529297 

Right  Ascen  at  Week(rad):  -0.2545632839E+001 

Argument  of  Perigee(rad):  -2.099498391 

Mean  Anom(rad):  -0.93458 18758E+000 

Af0(s):  -0.1583099365E-003 

Afl(s/s):  0.3637978807E-011 

week:  493 

********  Week  493  almanac  for  PRN-15  ******** 

ID:  15 

Health:  000 

Eccentricity:  0.1409053802E-002 

Time  of  Applicability(s):  405504.0000 

Orbital  Inclination(rad):  0.9571590424 

Rate  of  Right  Ascen(r/s):  -0.8178176358E-008 

SQRT(A)  (m  1/2):  5153.528320 

Right  Ascen  at  Week(rad):  -0.2594598174E+001 

Argument  of  Perigee(rad):  -0.741997719 

Mean  Anom(rad):  -0.42932701 1 1E+000 

Af0(s):  -0.2527236938E-003 

Afl(s/s):  -0.3637978807E-01 1 

week:  493 

********  Week  493  almanac  for  PRN-16  ******** 
ID:  16 

Health:  000 

Eccentricity:  0.5049705505E-002 

Time  of  Applicability(s):  405504.0000 
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Orbital  Inclination(rad):  0.9677886963 
Rate  of  Right  Ascen(r/s):  -0.7625203580E-008 
SQRT(A)  (m  1/2):  5153.715332 

Right  Ascen  at  Week(rad):  -0.4652485847E+000 
Argument  of  Perigee(rad):  -0.385978222 
Mean  Anom(rad):  0.9201985598E+000 

AfO(s):  0. 772476 1963E-004 

Afl(s/s):  -0.3637978807E-011 

week:  493 

********  Week  493  almanac  for  PRN-17  ******** 

ID:  17 

Health:  000 

Eccentricity:  0.4197120667E-002 

Time  of  Applicability(s):  405504.0000 

Orbital  Inclination(rad):  0.9598731995 

Rate  of  Right  Ascen(r/s):  -0.7799826562E-008 

SQRT(A)  (m  1/2):  5153.578613 

Right  Ascen  at  Week(rad):  0.5762445927E+000 

Argument  of  Perigee(rad):  -2.750026464 

Mean  Anom(rad):  -0.1657938957E+000 

Af0(s):  0.45  776367 19E-004 

Afl(s/s):  0.0000000000E+000 

week:  493 

********  Week  493  almanac  for  PRN-18  ******** 

ID:  18 

Health:  000 

Eccentricity:  0.9917736053E-002 

Time  of  Applicability(s):  405504.0000 

Orbital  Inclination(rad):  0.9440841675 

Rate  of  Right  Ascen(r/s):  -0.7879862096E-008 

SQRT(A)  (m  1/2):  5153.706543 

Right  Ascen  at  Week(rad):  0.2687054992E+001 

Argument  of  Perigee(rad):  -2.438257575 

Mean  Anom(rad):  0.1048157930E+001 
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AfO(s): 

Afl(s/s): 

week: 


-0.8583068848E-004 


0.3637978807E-011 
493 

********  Week  493  almanac  for  PRN-19  ******** 

ID:  19 

Health:  000 

Eccentricity:  0.5 12599945  IE-002 

Time  of  Applicability(s):  405504.0000 

Orbital  Inclination(rad):  0.9575653076 

Rate  of  Right  Ascen(r/s):  -0.7843482308E-008 

SQRT(A)  (m  1/2):  5153.721680 

Right  Ascen  at  Week(rad):  0.6326053 143 E+000 

Argument  of  Perigee(rad):  -0.391284585 

Mean  Anom(rad):  -0.9074705839E+000 

Af0(s):  0.3242492676E-004 

Afl(s/s):  0.0000000000E+000 

week:  493 

********  Week  493  almanac  for  PRN-20  ******** 

ID:  20 

Health:  000 

Eccentricity:  0.3940582275E-002 

Time  of  Applicability(s):  405504.0000 

Orbital  Inclination(rad):  0.9436454773 

Rate  of  Right  Ascen(r/s):  -0.7901689969E-008 

SQRT(A)  (m  1/2):  5153.634277 

Right  Ascen  at  Week(rad):  0.2633304238E+001 

Argument  of  Perigee(rad):  1.314607620 

Mean  Anom(rad):  0.1 195346236E+001 

Af0(s):  0.9155273438E-004 

Afl(s/s):  0.0000000000E+000 

week:  493 

********  Week  493  almanac  for  PRN-21  ******** 
ID:  21 


61 


Health:  000 

Eccentricity:  0.1462936401E-001 

Time  of  Applicability(s):  405504.0000 
Orbital  Inclination(rad):  0.9342021942 
Rate  of  Right  Ascen(r/s):  -0.8087226888E-008 
SQRT(A)  (m  1/2):  5153.579590 

Right  Ascen  at  Week(rad):  0.1628402948E+001 
Argument  of  Perigee(rad):  -2.648899078 
Mean  Anom(rad):  0.2324406743E+001 

Af0(s):  0.3051 7578 12E-004 

Afl(s/s):  0.0000000000E+000 

week:  493 

********  Week  493  almanac  for  PRN-22  ******** 

ID:  22 

Health:  000 

Eccentricity:  0.4887580872E-002 

Time  of  Applicability(s):  405504.0000 

Orbital  Inclination(rad):  0.9417476654 

Rate  of  Right  Ascen(r/s):  -0.7898051990E-008 

SQRT(A)  (m  1/2):  5153.617676 

Right  Ascen  at  Week(rad):  0.2692983747E+001 

Argument  of  Perigee(rad):  -1.823396802 

Mean  Anom(rad):  -0.7659566402E-001 

Af0(s):  0.2059936523E-003 

Afl(s/s):  0.0000000000E+000 

week:  493 

********  Week  493  almanac  for  PRN-23  ******** 

ID:  23 

Health:  000 

Eccentricity:  0.5777835846E-002 

Time  of  Applicability(s):  405504.0000 

Orbital  Inclination(rad):  0.9710483551 

Rate  of  Right  Ascen(r/s):  -0.8054485079E-008 

SQRT(A)  (m  1/2):  5153.600586 
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Right  Ascen  at  Week(rad):  -0.2576191902E+001 
Argument  of  Perigee(rad):  2.882402658 
Mean  Anom(rad):  -0.1446610928E+001 

AfO(s):  0.3  86238098  IE-003 

Afl(s/s):  0.0000000000E+000 

week:  493 

********  Week  493  almanac  for  PRN-24  ******** 

ID:  24 

Health:  000 

Eccentricity:  0.6956 100464E-002 

Time  of  Applicability(s):  405504.0000 

Orbital  Inclination(rad):  0.9508419037 

Rate  of  Right  Ascen(r/s):  -0.7930793799E-008 

SQRT(A)  (m  1/2):  5153.608398 

Right  Ascen  at  Week(rad):  0.1658670664E+001 

Argument  of  Perigee(rad):  -0.639523149 

Mean  Anom(rad):  0.9464823008E+000 

Af0(s):  0. 1 63078308  IE-003 

Afl(s/s):  0.3637978807E-01 1 

week:  493 

********  Week  493  almanac  for  PRN-25  ******** 

ID:  25 

Health:  000 

Eccentricity:  0.1179742813E-001 

Time  of  Applicability(s):  405504.0000 

Orbital  Inclination(rad):  0.9679870605 

Rate  of  Right  Ascen(r/s):  -0.8098 140825E-008 

SQRT(A)  (m  1/2):  5153.625488 

Right  Ascen  at  Week(rad):  -0.1628400803E+001 

Argument  of  Perigee(rad):  -1.22941 1483 

Mean  Anom(rad):  0. 146781361 1E+001 

Af0(s):  0.223 1597900E-003 

Afl(s/s):  0. 2546585 165E-0 10 

week:  493 
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********  Week  493  almanac  for  PRN-26  ******** 

ID:  26 

Health:  000 

Eccentricity:  0.1954650879E-001 

Time  of  Applicability(s):  405504.0000 

Orbital  Inclination(rad):  0.9932918549 

Rate  of  Right  Ascen(r/s):  -0.7832568372E-008 

SQRT(A)  (m  1/2):  5153.601074 

Right  Ascen  at  Week(rad):  -0.2528931379E+001 

Argument  of  Perigee(rad):  0.967672229 

Mean  Anom(rad):  -0.1885770559E+001 

Af0(s):  0.4 1 96166992E-004 

Afl(s/s):  0.3637978807E-011 

week:  493 

********  Week  493  almanac  for  PRN-27  ******** 

ID:  27 

Health:  000 

Eccentricity:  0.2107429504E-001 

Time  of  Applicability(s):  405504.0000 

Orbital  Inclination(rad):  0.9726352692 

Rate  of  Right  Ascen(r/s):  -0.797081 1566E-008 

SQRT(A)  (m  1/2):  5155.133789 

Right  Ascen  at  Week(rad):  -0.1592519999E+001 

Argument  of  Perigee(rad):  -1.651 172996 

Mean  Anom(rad):  0.5396271944E+000 

Af0(s):  0.2193450928E-004 

Afl(s/s):  0.3637978807E-01 1 

week:  493 

********  Week  493  almanac  for  PRN-28  ******** 

ID:  28 

Health:  000 

Eccentricity:  0. 1463699341E-001 

Time  of  Applicability(s):  405504.0000 
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Orbital  Inclination(rad):  0.9659786224 
Rate  of  Right  Ascen(r/s):  -0.7614289643E-008 
SQRT(A)  (m  1/2):  5153.596191 

Right  Ascen  at  Week(rad):  -0.4569123983E+000 
Argument  of  Perigee(rad):  -2.057136774 
Mean  Anom(rad):  0.2584432364E+000 

AfO(s):  -0.23  84 1 8579  IE-004 

Afl(s/s):  0.0000000000E+000 

week:  493 

********  Week  493  almanac  for  PRN-29  ******** 

ID:  29 

Health:  000 

Eccentricity:  0.3262042999E-002 

Time  of  Applicability(s):  405504.0000 

Orbital  Inclination(rad):  0.9600524902 

Rate  of  Right  Ascen(r/s):  -0.7814378478E-008 

SQRT(A)  (m  1/2):  5153.684082 

Right  Ascen  at  Week(rad):  0.5844268799E+000 

Argument  of  Perigee(rad):  -1.343433499 

Mean  Anom(rad):  0.2377104282E+001 

Af0(s):  -0.123977661  IE-004 

Afl(s/s):  0.3637978807E-011 

week:  493 

********  Week  493  almanac  for  PRN-30  ******** 

ID:  30 

Health:  000 

Eccentricity:  0.1 1 17467880E-001 

Time  of  Applicability(s):  405504.0000 

Orbital  Inclination(rad):  0.9512081146 

Rate  of  Right  Ascen(r/s):  -0.7788912626E-008 

SQRT(A)  (m  1/2):  5153.597656 

Right  Ascen  at  Week(rad):  -0.5279721022E+000 

Argument  of  Perigee(rad):  1.413174987 

Mean  Anom(rad):  0.6904801 130E+000 
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AfO(s): 

Afl(s/s): 

week: 


0. 12493 13354E-003 


0.3637978807E-01 1 
493 

********  Week  493  almanac  for  PRN-31  ******** 

ID:  31 

Health:  000 

Eccentricity:  0.7188796997E-002 

Time  of  Applicability(s):  405504.0000 

Orbital  Inclination(rad):  0.9710483551 

Rate  of  Right  Ascen(r/s):  -0.8021743270E-008 

SQRT(A)  (m  1/2):  5153.591309 

Right  Ascen  at  Week(rad):  -0.1536291 122E+001 

Argument  of  Perigee(rad):  -1.193828821 

Mean  Anom(rad):  -0.3108037591E+001 

Af0(s):  -0.5054473877E-004 

Afl(s/s):  0.0000000000E+000 

week:  493 

********  Week  493  almanac  for  PRN-32  ******** 

ID:  32 

Health:  000 

Eccentricity:  0.13553 14255E-001 

Time  of  Applicability(s):  405504.0000 

Orbital  Inclination(rad):  0.9657154083 

Rate  of  Right  Ascen(r/s):  -0.7672497304E-008 

SQRT(A)  (m  1/2):  5153.600586 

Right  Ascen  at  Week(rad):  0.2740463018E+001 

Argument  of  Perigee(rad):  -1.201977730 

Mean  Anom(rad):  -0.2252938390E+001 

Af0(s):  0.2975463867E-003 

Afl(s/s):  -0.3637978807E-011 

week:  493 
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APPENDIX  B.  VISUAL  C++  CODES  FOR  DOP  CALCULATION 


A.  DOPCALCULATORDLG.CPP 

// DOPCalculatorDlg.cpp  :  implementation  file 

#include  <tchar.h> 

#include  <urlmon.h> 

#pragma  comment(lib,  "urlmon.lib") 

#include  "stdafx.h" 

#include  <math.h> 

#include  <stdlib.h> 

#include  <string> 

#include  <iostream> 

#include  <fstream> 

#include  "DOP  Calculator.h" 

#include  "DOP  CalculatorDlg.h" 

#include  "Matrix. h" 

using  namespace  std; 

#ifdef  DEBUG 

#define  new  DEBUGNEW 

#endif 

//  CAboutDlg  dialog  used  for  App  About 

class  CAboutDlg  :  public  CDialog 

{ 

public: 

CAboutDlgQ; 

//  Dialog  Data 

enum  { IDD  =  IDD  ABOUTBOX  }; 
protected: 
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virtual  void  DoDataExchange(CDataExchange*  pDX);  //  DDX/DDV  support 


//  Implementation 
protected: 

DECLARE_MESSAGE_MAP() 

}; 


CAboutDlg: :CAboutDlg() :  CDialog(CAboutDlg::IDD) 

{ 

} 


void  CAboutDlg:  :DoDataExchange(CDataExchange*  pDX) 

{ 

CDialog::DoDataExchange(pDX); 

} 


BEGIN_MESSAGE_MAP(CAboutDlg,  CDialog) 

END_MESSAGE_MAP() 

//  CDOP  CalculatorDlg  dialog 

CDOP_CalculatorDlg::CDOP_CalculatorDlg(CWnd*  pParent  /*=NULL*/) 
:  CDialog(CDOP_CalculatorDlg::IDD,  pParent) 

,  yr(2009) 

,  mth(2) 

,  day(3) 

,  hr(0) 

,  mi(0) 

,  timezone(O) 

,  longi(O) 

,  lat(O) 

,  alt(O) 

,  obs(10) 

,  xdopO(_T("")) 

,  ydopO(_T("")) 

,  vdopO(_T("")) 


69 


,  tdopO(_T("")) 

,  hdopO(_T("")) 

,  pdopO(_T("")) 

,  gdopO(_T("")) 

,  daylightsaving(FALSE) 

{ 

mhlcon  =  AfxGetApp()->LoadIcon(IDR_MAINFRAME); 
m_pGraph=NULL; 


void  CDOP_CalculatorDlg::DoDataExchange(CDataExchange*  pDX) 

{ 

CDialog::DoDataExchange(pDX); 


DDX_Text(pDX,  IDC_YR,  yr); 
DDX_Text(pDX,  IDC_MTH,  mth); 
DDV_MinMaxDouble(pDX,  mth,  1,  12); 
DDX_Text(pDX,  IDCDAY,  day); 
DDV_MinMaxDouble(pDX,  day,  1,  31); 
DDX_Text(pDX,  IDC_HR,  hr); 
DDV_MinMaxDouble(pDX,  hr,  0,  23); 
DDX_Text(pDX,  IDC_MI,  mi); 
DDV_MinMaxDouble(pDX,  mi,  0,  59); 
DDX_Text(pDX,  IDC_TIMEZONE,  timezone); 
DDX_Text(pDX,  IDC_LONGI,  longi); 
DDV_MinMaxDouble(pDX,  longi,  -180,  180); 
DDX_Text(pDX,  IDC  LAT,  lat); 
DDV_MinMaxDouble(pDX,  lat,  -90,  90); 
DDX_Text(pDX,  IDC  ALT,  alt); 
DDX_Text(pDX,  IDC_OBS,  obs); 
DDX_Text(pDX,  IDC_XDOP0,  xdopO); 
DDX_Text(pDX,  IDC_YDOP0,  ydopO); 
DDX_Text(pDX,  IDC_VDOP0,  vdopO); 
DDX_Text(pDX,  IDC_TDOP0,  tdopO); 
DDX_Text(pDX,  IDC_HDOP0,  hdopO); 
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DDX_Text(pDX,  IDC_PDOPO,  pdopO); 

DDX_Text(pDX,  IDC_GDOPO,  gdopO); 

DDX_Check(pDX,  IDC  D A YLIGHT SAY,  daylightsaving); 


BEGIN_MESSAGE_MAP(CDOP_CalculatorDlg,  CDialog) 

ONWMSYSCOMMANDO 
ONWMP  AINT  () 

ONWMQUERYDRAGICONO 
//} }  AFX_MSG_MAP 

ON_BN_CLICKED(IDC_ALMANAC,  &CDOP_CalculatorDlg::OnBnClickedAlmanac) 
ON_BN_CLICKED(IDC_COMPUTE,  &CDOP_CalculatorDlg::OnBnClickedCompute) 
ON_BN_CLICKED(IDC_XDOP_RADIO,  &CDOP_CalculatorDlg::OnBnClickedXdopRadio) 
ON_BN_CLICKED(IDC_YDOP_RADIO,  &CDOP_CalculatorDlg::OnBnClickedYdopRadio) 
ON_BN_CLICKED(IDC_VDOP_RADIO,  &CDOP_CalculatorDlg::OnBnClickedVdopRadio) 
ON_BN_CLICKED(IDC_TDOP_RADIO,  &CDOP_CalculatorDlg::OnBnClickedTdopRadio) 
ON_BN_CLICKED(IDC_HDOP_RADIO,  &CDOP_CalculatorDlg::OnBnClickedHdopRadio) 
ON_BN_CLICKED(IDC_PDOP_RADIO,  &CDOP_CalculatorDlg::OnBnClickedPdopRadio) 
ON_BN_CLICKED(IDC_GDOP_RADIO,  &CDOP_CalculatorDlg::OnBnClickedGdopRadio) 
END_MESSAGE_MAP() 

//  CDOP  CalculatorDlg  message  handlers 

BOOL  CDOP  CalculatorDlg:  :OnInitDialog() 

{ 

CDialog::OnInitDialog(); 

//  Add  "About..."  menu  item  to  system  menu. 

//  IDM  ABOUTBOX  must  be  in  the  system  command  range. 

ASSERT((IDM_ABOUTBOX  &  OxFFFO)  ==  IDM  ABOUTBOX); 
ASSERT(IDM_ABOUTBOX  <  OxFOOO); 

CMenu*  pSysMenu  =  GetSystemMenu(FALSE); 
if  (pSysMenu  !=NULL) 

{ 
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CString  strAboutMenu; 

strAboutMenu. LoadString(IDSABOUTBOX); 
if  (!  strAboutMenu.  IsEmptyi)) 

{ 

pSysMenu->AppendMenu(MF_SEPARATOR); 

pSysMenu->AppendMenu(MF_STRING,  IDMABOUTBOX,  strAboutMenu); 

} 


//  Set  the  icon  for  this  dialog.  The  framework  does  this  automatically 
//  when  the  application's  main  window  is  not  a  dialog 
SetIcon(m_hIcon,  TRUE);  //  Set  big  icon 

SetIcon(m_hIcon,  FALSE);  //  Set  small  icon 

//  TODO:  Add  extra  initialization  here 
//This  change  the  position  of  default  graph 

m_pGraph=newCGraph(this,230,230,300,250,G_WHITESCHEME); 

mlpfs^NULL; 

m_pPlotItems=NULL; 

almanaclocation  =  "FPWMyDocsWVisual  Studio  2005WProjects\\Test3\\Test3\\current.al3"; 
gotdata  =  false; 

return  TRUE;  //  return  TRUE  unless  you  set  the  focus  to  a  control 

} 

void  CDOP_CalculatorDlg::OnSysCommand(UINT  nID,  LPARAM  IParam) 

{ 

if  ((nID  &  OxFFFO)  ==  IDM  ABOUTBOX) 

{ 

CAboutDlg  dlgAbout; 
dlgAbout.DoModalQ; 

} 

else 

{ 

CDialog::OnSysCommand(nID,  IParam); 


72 


//  If  you  add  a  minimize  button  to  your  dialog,  you  will  need  the  code  below 
//  to  draw  the  icon.  For  MFC  applications  using  the  document/view  model, 
//  this  is  automatically  done  for  you  by  the  framework. 


void  CDOP  CalculatorDlg:  :OnPaint() 

{ 

if  (IsIconicO) 

{ 

CPaintDC  dc(this);  //  device  context  for  painting 

S  endMessage(  WMICONERA  SEBKGND , 
reinterpret_cast<WPARAM>(dc.GetSafeHdc()),  0); 

//  Center  icon  in  client  rectangle 
int  cxlcon  =  GetSystemMetrics(SMCXICON); 
int  cylcon  =  GetSystemMetrics(SMCYICON); 
CRect  rect; 

GetClientRect(&rect); 

int  x  =  (rect.Width()  -  cxlcon  +  1)  /  2; 

int  y  =  (rect.Height()  -  cylcon  +  1)  /  2; 

//  Draw  the  icon 
dc.DrawIcon(x,  y,  m  hlcon); 

} 

else 

{ 

CDialog:  :OnPaint(); 
m_pGraph->PaintGraph( ) ; 

} 

} 


//  The  system  calls  this  function  to  obtain  the  cursor  to  display  while  the  user  drags 
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//  the  minimized  window. 

HCURSOR  CDOPCalculatorDlg:  :OnQueryDragIcon() 

{ 

return  static_cast<HCURSOR>(m_hIcon); 

} 


void  CDOP_CalculatorDlg::OnBnClickedAlmanac() 

{ 

//Downloads  latest  almanac  data  from  the  USCG  website  to  specified  location 

if(URLDownloadToFile(NULL,  _T("http://www. navcen.uscg.gov/gps/ciirrent/ciirrent.al3"), 
almanaclocation,  0,  NULL  )) 

{ 

AfxMessageBox((CString)"Error:  File  Not  Found  or  you  are 
offline", MBICONERROR); 

} 

else 

{ 

AfxMessageBox((CString)"Latest  Almanac  Downloaded", MB  OK); 

} 

} 


void  CDOP  CalculatorDlg:  :OnBnClickedCompute() 

{ 

//Clears  any  existing  plot  data 

m_lpfs=NULL; 

m_pPlotItems=NULL; 

//Updates  user  input 
UpdateData(TRUE); 

if  (daylightsaving  ==  true) 
daylightsav  =  1 ; 

else 

daylightsav  =  0; 
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// - 

//Convert  Target's  Latitude  ,  Longitude  and  Altitude  to  ECEF  Coordinates 

// - 

//WGS84  ellopsoid  constants 

double  a  =  6378137; 

double  es  =  8.1819190842622e-2; 

double  sealevel  =  0;  //Altitude  of  sealevel  is  zero  meter 

double  pi  =  3.1415926535898; 
double  lat  rad  =  lat  *  pi/180; 
double  long_rad  =  longi  *  pi/180; 

double  N  =  a/sqrt(  1  -  (  pow(es,  2) )  *  (  pow(  sin(lat_rad),  2  ) ) );  //Prime  vertical  radius  of 
curvature 

double  xTgt  =  (N  +  alt)  *  cos(lat  rad)  *  cos(long  rad);  //Target  x  ECEF  coordinates  (meters) 
double  yTgt  =  (N  +  alt)  *  cos(lat_rad)  *  sin(long_rad);  //Target  y  ECEF  coordinates  (meters) 
double  zTgt  =  ( ( (1  -  pow(es,  2) )  *  N)  +  alt )  *  sin(lat_rad);  //Target  z  ECEF  coordinates  (meters) 

double  xLocalRef  =  (N  +  sealevel)  *  cos(lat  rad)  *  cos(long  rad);  //x  ECEF  coordinates  for  Local 
Reference  Pt  for  ENU  frame  (meters) 

double  yLocalRef  =  (N  +  sealevel)  *  cos(lat  rad)  *  sin(long  rad);  //y  ECEF  coordinates  for  Local 
Reference  Pt  for  ENU  frame  (meters) 

double  zLocalRef  =  ( ( (1  -  pow(es,  2) )  *  N)  +  sealevel )  *  sin(lat_rad);  //z  ECEF  coordinates  for 
Local  Reference  Pt  for  ENU  frame  (meters) 

// - 

//Convert  Almanac  Data  to  Satellite  Position  in  East-North-Up  (ENU)  Coordinates  and  Check  the 
Line  of  Sight 

// - 

//Constants 

double  io  =  0.3  *  pi;  //Inclination  angle  @  ref.  time  (rad) 

double  mju  =  3.986005el4;  //WGS84  value  of  the  Earth's  universal  gravitational  parameter  for 
GPS  user  (meterA3/secA2) 

double  OMEGAdote  =  7.2921 151467e-5;  //WGS84  value  of  the  Earth's  rotation  rate  (rad/sec) 
double  c  =  2.99792458e8;  //Speed  of  light  (meter/sec) 

//Reads  almanac  data  from  Almanac  file 
char  inputline[100]; 
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inputline[0]='\0'; 


int  NumSV ;  //Number  of  satellites 
int  wn;  //GPS  week  no. 

int  toa;  //Time  of  applicability  of  Almanac  (sec)  (range:  0  to  604,784) 

//Reads  the  almanac  data  from  specified  location 

ifstream  in(almanac_location,  ios::in);  //opens  the  almanac  file 

in  »  NumSV  »  inputline  »  wn  »  toa;  //reads  time  info  from  almanac  file 

//Convert  date  to  GPS  time 

double  JAN61980  =  44244  +  2400000.5;  //Date  when  GPS  started 
double  sec_per_day  =  86400; 

double  JD  =  367*yr  -  floor(7*(yr  +  floor((mth  +  9)/12))/4)  -  floor(3*(floor((yr  +  (mth-9)/7)/100)  + 
l)/4)  +  floor(275*mth/9)  +  day  +  1721028.5; 

double  fJD  =  (mi/1440)  +  (hr/24); 

JD  =  JD  +  fJD  -  (timezone  +  daylightsav)/24; 

//Correct  for  leap  seconds  up  till  Mar  2009 
int  GpsLeapSec; 

if  ((JD>=2444786.5)&&(JD<2445 151.5)) 

GpsLeapSec=l; 

else  if  ((JD>=2445 1 5 1 .5)&&(JD<24455 16.5)) 

GpsLeapSec=2; 

else  if  ((JD>=24455 16.5)&&(JD<2446247.5)) 

GpsLeapSec=3; 

else  if  ((JD>=2446247.5)&&(JD<2447161 .5)) 

GpsLeapSec=4; 

else  if((JD>=2447161.5)&&(JD<2447892.5)) 

GpsLeapSec=5; 

else  if  ((JD>=2447892.5)&&(JD<2448257.5)) 

GpsLeapSec=6; 

else  if  ((JD>=2448257.5)&&(JD<2448804.5)) 

GpsLeapSec=7; 
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else  if((JD>=2448804.5)&&(JD<2449169.5)) 

GpsLeapSec=8; 

else  if((JD>=2449169.5)&&(JD<2449534.5)) 

GpsLeapSec=9; 

else  if((JD>=2449534.5)&&(JD<2450083.5)) 

GpsLeapSec=10; 

else  if((JD>=2450083.5)&&(JD<2450630.5)) 

GpsLeapSec=ll; 

else  if  ((JD>=2450630.5)&&(JD<245 1 179.5)) 

GpsLeapSec=12; 

else  if  ((JD>=2451 179.5)&&(JD<2453371.5)) 

GpsLeapSec=13; 

else  if  ((JD>=2453371.5)&&(JD<2454466.5)) 

GpsLeapSec=14; 

else  GpsLeapSec=15;  //Add  when  new  seconds  are  added 

JD  =  JD  +  GpsLeapSec/sec_per_day; 

double  gps_week  =  floor/  (JD  -  JAN61980)/7  ); 

double  sec_of_week  =  ( (JD  -  JAN61980)  -  gps_week  *  7  )  *  sec_per_day; 

//End  of  GPS  time  conversion 

while  (gps_week  >=  1024) 

{ 

gps_week  =  gps_week  -  1024; 

} 


double  tk  =  (gps_week  -  wn)  *  604800  +  (sec_of_week  -  toa);  //Time  since  toa  (sec)  (range 
302400  to  302400) 

//Need  to  increase  the  array  size  if  there  are  more  than  3 1  GPS  satellites 

int  PRN; 

int  SVN; 

int  URA; 

double  ec; 

double  del  ik; 

double  OMEGAdot; 
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double  sqrtA; 

double  OMEGAo; 

double  omega; 

double  Mo; 

double  AfO; 

double  Afl ; 

double  SVHEALTH; 

double  SVCONFIG; 

double  EAST[31],  NORTH[31],  UP[31]; 

int  LOS[31]; 

int  il; 

for  (i  1=0;  il<NumSV;  il++)  //reads  data  for  different  satellites 

{ 

//Read  the  data  from  Almanac  file 
in  »  PRN  »  SVN  »  URA; 
in  »  ec  »  delik  »  OMEGAdot; 
in  »  sqrtA  »  OMEGAo  »  omega; 
in  »  Mo  »  AfD  »  Afl ; 
in  »  SVHEALTH  »  SVCONFIG; 

//Conversion  of  some  almanac  data  to  radians 
del  ik  =  del  ik  *  pi; 

OMEGAdot  =  OMEGAdot  *  pi; 

OMEGAo  =  OMEGAo  *  pi; 
omega  =  omega  *  pi; 

Mo  =  Mo  *  pi; 

//Position  of  Satellite  Calculations 

double  A  =  pow(sqrtA,  2); 

double  n  =  sqrt(mju  /  pow(A,  3)); 

double  Mk  =  Mo  +  (tk  *  n);  //Mean  anomaly  (rad) 

//Start  values  for  iterative  solution  of  Kepler  Equation 
double  Ek  =  Mk; 
double  Eold  =  0; 
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while  (  abs(Ek  -  Eold)  >=  le-10  ) 

{ 

Eold  =  Ek; 

Ek  =  Mk  +  ec*sin(Ek); 

} 

double  vk  =  atan2(  (  sqrt(  1  -  pow(ec,  2) )  *  sin(Ek) )/(  1  -  ec  *  cos(Ek) ),  (  cos(Ek)  -  ec  )/( 
1  -  ec  *  cos(Ek)  ) );  //True  Anomaly  (rad) 

Ek  =  acos(  (ec  +  cos(vk) )  /  (  1  +  ec  *  cos(vk) ) ); 

double  uk  =  omega  +  vk;  //Argument  of  latitude  (rad) 

double  rk  =  A  *  (  1  -  ec  *  cos(Ek) );  //Corrected  radius  (meter) 

double  ik  =  io  +  del  ik;  //Corrected  inclination  (rad) 

double  xkl  =  rk  *  cos(uk);  //x  position  in  orital  plane  (meter) 

double  ykl  =  rk  *  sin(uk);  l/y  position  in  orital  plane  (meter) 

double  OMEGAk  =  OMEGAo  +  ( (OMEGAdot  -  OMEGAdote)  *  tk  )  -  (OMEGAdote  * 
toa);  //Corrected  longitude  of  ascending  node  (rad) 

//Calculations  for  ECEF  coordinates 

double  xk  =  (  xkl  *  cos(OMEGAk)  )  -  (  ykl  *  cos(ik)  *  sin(OMEGAk) );  //Satellite  x 
ECEF  coordinate  (meter) 

double  yk  =  (  xkl  *  sin(OMEGAk) )  +  ( ykl  *  cos(ik)  *  cos(OMEGAk) );  //Satellite  y 
ECEF  coordinate  (meter) 

double  zk  =  ykl  *  sin(ik);  //Satellite  z  ECEF  coordinate  (meter) 

//Convert  ECEF  coordinates  to  East-North-Up  (ENU)  coordinates 

double  East  =  -sin(long  rad)  *  (xk  -  xLocalRef)  +  cos(long  rad)  *  (yk  -  yLocalRef); 

double  North  =  ( -sin(lat  rad)  *  cos(long  rad)  *  (xk  -  xLocalRef)  )  -  (  sin(lat  rad)  * 
sin(long  rad)  *  (yk  -  yLocalRef)  )  +  (  cos(lat  rad)  *  (zk  -  zLocalRef)  ); 

double  Up  =  (  cos(lat  rad)  *  cos(long  rad)  *  (xk  -  xLocalRef)  )  +  (  cos(lat  rad)  * 
sin(long  rad)  *  (yk  -  yLocalRef)  )  +  (  sin(lat  rad)  *  (zk  -  zLocalRef)  ); 

EAST[il]  =  East; 

NORTH[il]  =  North; 

UP[il]  =  Up; 

//End  of  Position  Calculations  of  Satellites 


//Determine  Line  of  Sight  between  Target  and  Satellite 


79 


double  mag_Tgt  =  sqrt(  pow(xTgt,  2)  +  pow(yTgt,  2)  +  pow(zTgt,  2) );  //Distance  of 
Target  from  Earth  center  (meter) 

double  xTgttoSV  =  xk  -  xTgt; 

double  yTgttoSV  =  yk  -  yTgt; 

double  zTgttoSV  =  zk  -  zTgt; 

double  mag_TgttoSV  =  sqrt(  pow(  (xk  -  xTgt),  2  )  +  pow(  (yk  -  yTgt),  2  )  +  pow(  (zk  - 
zTgt),  2  ) );  //Distance  from  Target  to  Satellite 

double  AngleFromTgt  =  acos(((xTgttoSV*xTgt)  +  (yTgttoSV*yTgt)  + 
(zTgttoSV*zTgt))/(mag_TgttoSV*mag_Tgt)); 

int  Los; 

if  (AngleFromTgt  <  (pi/2-(obs*pi/180))  &&  SVHEALTH  ==  0) 

Los  =  1;  //There  is  Line  of  Sight 

else 

Los  =  0;  //No  Line  of  Sight 
LOS[il]  =  Los; 

//End  of  Determination  of  Line  of  Sight 

} 


//Assign  Visible  Satellites 
int  i2; 
int  i3  =  0; 
double  SV[31][3]; 

for(i2=0;  i2<NumSV;  i2++) 

{ 

if  (LOS[i2]  ==  1) 

{ 

//Assigning  coordinates  to  respective  visible  satellites 
SV[i3][0]  =  EAST[i2]; 

SV[i3][l]  =  NORTFl[i2]; 

SV[i3][2]  =  UP[i2]; 
i3++; 


80 


int  NumVisibleSV  =  i3; 


// - 

//Calculate  DOP  for  Specified  Time 

// - 

//Pseudo-Range  and  Directional  Derivative  Loop 
int  i4; 

double  r[3 1],  Dx[3 1],  Dy[3 1],  Dz[3 1],  Dt[3 1]; 

for(i4=0;  i4<NumVisibleSV;  i4++) 

{ 

//Calculate  pseudo-ranges  from  target  position  to  visible  satellites 

r[i4]  =  sqrt(  pow(SV[i4][0],  2)  +  pow(SV[i4][l],  2)  +  pow((SV[i4][2]  -  alt),  2) ); 

//Calculate  directional  derivatives  for  East,  North,  Up  and  Time 

Dx[i4]  =  SV[i4][0]/r[i4]; 

Dy[i4]  =  SV[i4][l]/r[i4]; 

Dz[i4]  =  (S V[i4] [2]  -  alt)  /  r[i4]; 

Dt[i4]  =  -1; 

} 

//Produce  the  Covariance  Matrix  from  the  Directional  Derivatives 
int  i5,  i6; 

double  Alp[31][4]; 


for  (i5=0;  i5<3 1 ;  i5++) 

{ 


for  (i6=0;  i6<4;  i6++) 

{ 


Alp[i5][i6]  =  0;  //Initialize  Alp 


for  (i5=0;  i5<NumVisibleSV;  i5++) 

{ 

Alp[i5][0]  =  Dx[i5]; 
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Alp[i5][l]  =  Dy[i5]; 
Alp[i5][2]  =  Dz[i5]; 
Alp[i5][3]  =  Dt[i5]; 


//Transpose  Alp  to  get  Brv 
int  i7,  i8; 

double  Brv[4][31]; 


for  (i7=0;  i7<4;  i7++) 

{ 


for  (i8=0;  i8<3 1 ;  i8++) 

{ 


Alp[i5][i6]  =  0;  //Initialize  Brv 


} 

for  (i7=0;  i7<4;  i7++) 

{ 

i8  =  0; 

while  (i8<NumVisibleSV) 

{ 

Brv[i7][i8]  =  Alp[i8][i7]; 
i8++; 


//Matrix  multiplication  of  Brv  and  Alp 
int  i9,  ilO,  ill; 
double  Chl[4] [4] ; 

for  (i9=0;  i9<4;  i9++) 

{ 

for  (i  10=0;  il0<4;  il0++) 

{ 

Chl[i9][il0]  =  0;  //Initialize  Chi 
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} 


for  (i9=0;  i9<4;  i9++) 

{ 

for  (i  10=0;  il0<4;  il0++) 

{ 

for  (i  1 1=0;  ill<31;  ill++) 

Chl[i9][il0]  =  Chl[i9][il0]  +  Brv[i9][il  l]*Alp[il  l][il0]; 


//Inverse  Chi 
Matrix  m(4,4); 

int  i  12,  i  13; 

for  (i  12=0;  il2<4;  il2++) 

{ 

for  (i  1 3=0;  il3<4;  il3++) 

{ 

m(il2+l,  i  13+1)  =  Chl[il2][i  13];  //Assign  Chi  to  matrix,  m 


Matrix  minv(4,4); 

double  det=Matrix::inv(m,minv); 

double  Dlt[4] [4]; 

for  (i  12=0;  il2<4;  il2++) 

{ 

for  (i  13=0;  il3<4;  il3++) 

{ 

Dlt[il2][il3]  =  minv(il2+l,  il3+l);  //Assign  inversed  matrix,  minv,  to  Dlt 
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//Calculate  DOP 

xdop  =  sqrt(Dlt[0][0]); 

ydop  =  sqrt(Dlt[l][l]); 

vdop  =  sqrt(Dlt[2][2]); 

tdop  =  sqrt(Dlt[3][3]); 

hdop  =  sqrt(Dlt[0][0]  +  Dlt[l][l]); 

pdop  =  sqrt(Dlt[0][0]  +  Dlt[l][l]  +  Dlt[2][2]); 

gdop  =  sqrt(Dlt[0][0]  +  Dlt[l][l]  +  Dlt[2][2]  +  Dlt[3][3]); 

//Display  DOP  to  dialogue  box 
xdopO.Format((CString)"%.3f',  xdop); 
ydopO.Format((CString)"%.3f\  ydop); 
vdopO.Format((CString)"%.3f\  vdop); 
tdopO.Format((CString)"%.3f',  tdop); 
hdopO.Format((CString)"%.3f\  hdop); 
pdopO.Format((CString)"%.3f',  pdop); 
gdopO.Format((CString)"%.3f\  gdop); 

UpdateData(FALSE); 

//End  of  DOP  Calculation  for  Specified  Time 

// - 

//Calculate  DOP  for  24hr  Period  (hourly  sampling) 

// - 

double  xdop2,  ydop2,  zdop2,  tdop2,  gdop2,  pdop2,  hdop2;//  DOP  variables  for  graph  plotting 

//Check  if  there  is  any  existing  graph.  If  so,  delete  it. 
if(m_pPlotItems  !=NULL) 

{ 

delete  []m_pPlotItems; 
m_pPlotItems=NULL ; 

} 

if(m_lpfs  !=NULL) 

{ 

m_pGraph->ClearFunction( ) ; 
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delete  m  lpfs; 


} 

//Loop  to  calculate  DOP  values  every  lhr 
double  hrl  =  0; 
double  mil  =0; 

num_items  =  24;  //number  of  plot  points  on  the  x-axis 

XDOP_PlotPoints=new  double[num_items*2]; 

YDOP_PlotPoints=new  double[num_items*2]; 

VDOP_PlotPoints=new  double[num_items*2]; 

TDOP_PlotPoints=new  double[num_items*2]; 

HDOP_PlotPoints=new  double[num_items*2]; 

PDOP_PlotPoints=new  double[num_items*2]; 

GDOP_PlotPoints=new  double[num_items*2]; 

//Loop  for  hourly  sampling 

for(int  i  14=0;  il4<=num_items*2;  il4+=2) 

{ 

JD  =  367*yr  -  floor(7*(yr  +  floor((mth  +  9)/12))/4)  -  floor(3*(floor((yr  +  (mth-9)/7)/100) 
+  l)/4)  +  floor(275*mth/9)  +  day  +  1721028.5; 

double  fJD  =  (mil/1440)  +  (hr  1/24); 

JD  =  JD  +  fJD  -  (timezone  +  daylightsav)/24; 

//Correct  for  leap  seconds  up  till  Mar  2009 
int  GpsLeapSec; 

if  ((JD>=2444786.5)&&(JD<2445 151.5)) 

GpsLeapSec=l; 

else  if  ((JD>=2445 151 ,5)&&(JD<24455 1 6.5)) 

GpsLeapSec=2; 

else  if((JD>=2445516.5)&&(JD<2446247.5)) 

GpsLeapSec=3; 

else  if((JD>=2446247.5)&&(JD<2447161.5)) 


85 


GpsLeapSec=4; 

else  if((JD>=2447161.5)&&(JD<2447892.5)) 

GpsLeapSec=5; 

else  if  ((JD>=2447892.5)&&(JD<2448257.5)) 

GpsLeapSec=6; 

else  if  ((JD>=2448257.5)&&(JD<2448804.5)) 

GpsLeapSec=7; 

else  if((JD>=2448804.5)&&(JD<2449169.5)) 

GpsLeapSec=8; 

else  if((JD>=2449169.5)&&(JD<2449534.5)) 

GpsLeapSec=9; 

else  if((JD>=2449534.5)&&(JD<2450083.5)) 

GpsLeapSec=10; 

else  if((JD>=2450083.5)&&(JD<2450630.5)) 

GpsLeapSec=ll; 

else  if  ((JD>=2450630.5)&&(JD<245 1 179.5)) 

GpsLeapSec=12; 

else  if  ((JD>=2451 179.5)&&(JD<2453371.5)) 

GpsLeapSec=13; 

else  if((JD>=2453371.5)&&(JD<2454466.5)) 

GpsLeapSec=14; 

else  GpsLeapSec=15;  //Add  when  new  seconds  are  added 

JD  =  JD  +  GpsLeapSec/sec_per_day; 
double  gps_week  =  floor/  (JD  -  JAN61980)/7  ); 

double  sec_of_week  =  ( (JD  -  JAN61980)  -  gps_week  *  7  )  *  sec_per_day; 

// End  of  GPS  time  conversion 

while  (gps_week  >=  1024) 

{ 

gps_week  =  gps_week  -  1024; 

} 

//Reads  the  almanac  data  from  specified  location 

ifstream  in(almanac_location,  ios::in);  //opens  the  almanac  file 

in  »  NumSV  »  inputline  »  wn  »  toa;  //reads  time  info  from  almanac  file 
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double  tk  =  (gps_week  -  wn)  *  604800  +  (sec_of_week  -  toa);  //Time  since  toa  (sec) 
(range:  -302400  to  302400) 

//Need  to  increase  the  array  size  if  there  are  more  than  31  GPS  satellites 

int  PRN; 

int  SVN; 

int  URA; 

double  ec; 

double  delik; 

double  OMEGAdot; 

double  sqrtA; 

double  OMEGAo; 

double  omega; 

double  Mo; 

double  AfO; 

double  Afl ; 

double  SVHEALTH; 

double  SVCONFIG; 

double  EAST[3 1],  NORTH[3 1],  UP[3 1]; 

int  LOS[3 1]; 

int  il; 

for  (il=0;  il<NumSV;  il++)  //reads  data  for  different  satellites 

{ 

//Read  the  data  from  Almanac  file 
in  »  PRN  »  SVN  »  URA; 
in  »  ec  »  del  ik  »  OMEGAdot; 
in  »  sqrtA  »  OMEGAo  »  omega; 
in  »  Mo  »  AfO  »  Afl ; 
in  »  SVHEALTH  »  SVCONFIG; 

//Conversion  of  some  almanac  data  to  radians 
del  ik  =  del  ik  *  pi; 

OMEGAdot  =  OMEGAdot  *  pi; 

OMEGAo  =  OMEGAo  *  pi; 
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omega  =  omega  *  pi; 

Mo  =  Mo  *  pi; 

//Position  of  Satellite  Calculations 

double  A  =  pow(sqrtA,  2); 

double  n  =  sqrt(mju  /  pow(A,  3)); 

double  Mk  =  Mo  +  (tk  *  n);  //Mean  anomaly  (rad) 

//Start  values  for  iterative  solution  of  Kepler  Equation 
double  Ek  =  Mk; 
double  Eold  =  0; 

while  (  abs(Ek  -  Eold)  >=  le-10  ) 

{ 

Eold  =  Ek; 

Ek  =  Mk  +  ec*sin(Ek); 

} 


double  vk  =  atan2(  (  sqrt(  1  -  pow(ec,  2) )  *  sin(Ek) )/(  1  -  ec  *  cos(Ek) ),  ( 
cos(Ek)  -  ec  )/(  1  -  ec  *  cos(Ek) ) );  //True  Anomaly  (rad) 

Ek  =  acos(  (ec  +  cos(vk) )  /  (  1  +  ec  *  cos(vk) ) ); 

double  uk  =  omega  +  vk;  //Argument  of  latitude  (rad) 

double  rk  =  A  *  (  1  -  ec  *  cos(Ek) );  //Corrected  radius  (meter) 

double  ik  =  io  +  del  ik;  //Corrected  inclination  (rad) 

double  xkl  =  rk  *  cos(uk);  //x  position  in  orital  plane  (meter) 

double  ykl  =  rk  *  sin(uk);  l/y  position  in  orital  plane  (meter) 

double  OMEGAk  =  OMEGAo  +  ( (OMEGAdot  -  OMEGAdote)  *  tk )  - 
(OMEGAdote  *  toa);  //Corrected  longitude  of  ascending  node  (rad) 

//Calculations  for  ECEF  coordinates 

double  xk  =  (  xkl  *  cos(OMEGAk) )  -  (  ykl  *  cos(ik)  *  sin(OMEGAk) ); 
//Satellite  x  ECEF  coordinate  (meter) 

double  yk  =  (  xkl  *  sin(OMEGAk) )  +  (  ykl  *  cos(ik)  *  cos(OMEGAk) ); 
//Satellite  y  ECEF  coordinate  (meter) 

double  zk  =  ykl  *  sin(ik);  //Satellite  z  ECEF  coordinate  (meter) 

//Convert  ECEF  coordinates  to  East-North-Up  (ENU)  coordinates 
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yLocalRef); 


double  East  =  -sin(longrad)  *  (xk  -  xLocalRef)  +  cos(long  rad)  *  (yk  - 


double  North  =  ( -sin(latrad)  *  cos(longrad)  *  (xk  -  xLocalRef) )  -  ( 
sin(latrad)  *  sin(long  rad)  *  (yk  -  yLocalRef)  )  +  (  cos(latrad)  *  (zk  -  zLocalRef)  ); 

double  Up  =  (  cos(latrad)  *  cos(longrad)  *  (xk  -  xLocalRef)  )  +  (  cos(latrad) 
*  sin(long  rad)  *  (yk  -  yLocalRef)  )  +  (  sin(latrad)  *  (zk  -  zLocalRef)  ); 


EAST[il]  =  East; 

NORTH[il]  =  North; 

UP[il]  =  Up; 

//End  of  Position  of  Satellite  Calculations 


//Determine  Line  of  Sight  between  Target  and  Satellite 

double  mag_Tgt  =  sqrt(  pow(xTgt,  2)  +  pow(yTgt,  2)  +  pow(zTgt,  2) ); 
//Distance  of  Target  from  Earth  center  (meter) 

double  mag_SV  =  sqrt(  pow(xk,  2)  +  pow(yk,  2)  +  pow(zk,  2) );  //Distance  of 
Satellite  from  Earth  center  (meter) 

double  AngleTOS  =  acos(  ( (xk  *  xTgt)  +  (yk  *  yTgt)  +  (zk  *  zTgt) )  /  (mag_Tgt 
*  mag_SV) );  //Angle  between  Target  and  Satellite  with  origin  at  center  of  Earth  (rad) 

double  mag_SVproj  =  mag_SV  *  cos(  AngleTOS);  //Magnitude  of  projection  of 
the  Satellite  vector  onto  Target  vector  (meter) 

double  xTgttoSV  =  xk  -  xTgt; 

double  yTgttoSV  =  yk  -  yTgt; 

double  zTgttoSV  =  zk  -  zTgt; 

double  xTgttoSVproj  =  ((mag_SVproj  -  mag_Tgt)  /  mag_Tgt)  *  xTgt; 

double  yTgttoSVproj  =  ((mag_SVproj  -  mag_Tgt)  /  mag_Tgt)  *  yTgt; 

double  zTgttoSVproj  =  ((mag  SVproj  -  magTgt)  /  mag  Tgt)  *  zTgt; 

double  mag_TgttoSV  =  sqrt(  pow(  (xk  -  xTgt),  2  )  +  pow(  (yk  -  yTgt),  2  )  + 
pow(  (zk  -  zTgt),  2  ) );  //Distance  from  Target  to  Satellite 

double  AngleFromTgt  =  acos(  ( (xTgttoSV  *  xTgttoSVproj)  +  (yTgttoSV  * 
yTgttoSVproj)  +  (zTgttoSV  *  zTgttoSVproj) )  /  ( (mag_SVproj  -  mag_Tgt)  *  mag_TgttoSV  ) );  //Angle 
between  Target  and  Satellite  with  origin  at  Target  (rad) 

int  Los; 

if  (  mag  SVproj  >  mag  Tgt  &&  AngleTOS  <  (pi/2)  &&  AngleFromTgt  <  (  pi/2 

-  (obs  *  pi/180) ) ) 

Los  =  1;  //There  is  Line  of  Sight 

else 

Los  =  0;  //No  Line  of  Sight 
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LOS[il]  =  Los; 

//End  of  Determination  of  Line  of  Sight 


//Assign  Visible  Satellites 
int  i2; 
int  i3  =  0; 
double  SV[31][3]; 

for(i2=0;  i2<NumSV;  i2++) 

{ 

if  (LOS[i2]  ==  1) 

{ 

//Assigning  coordinates  to  respective  visible  satellites 
SV[i3][0]  =  EAST[i2]; 

SV[i3][l]  =NORTH[i2]; 

SV[i3][2]  =  UP[i2]; 
i3++; 

} 

} 


int  NumVisibleSV  =  i3; 

// - 

//Calculate  DOP 

// 

//Pseudo-Range  and  Directional  Derivative  Loop 
int  i4; 

double  r[31],  Dx[31],  Dy[31],  Dz[31],  Dt[31]; 

for(i4=0;  i4<NumVisibleSV;  i4++) 

{ 

//Calculate  pseudo-ranges  from  target  position  to  visible  satellites 

r[i4]  =  sqrt(  pow(SV[i4][0],  2)  +  pow(SV[i4][l],  2)  +  pow((SV[i4][2]  -  alt),  2) ); 
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//Calculate  directional  derivatives  for  East,  North,  Up  and  Time 
Dx[i4]  =  SV[i4][0]/r[i4]; 

Dy[i4]  =  SV[i4][l]  /  r[i4] ; 

Dz[i4]  =  (S V[i4] [2]  -  alt)  /  r[i4]; 

Dt[i4]  =  -1; 

} 

//Produce  the  Covariance  Matrix  from  the  Directional  Derivatives 
int  i5,  i6; 

double  Alp[31][4]; 


for  (i5=0;  i5<3 1 ;  i5++) 

{ 


for  (i6=0;  i6<4;  i6++) 

{ 


Alp[i5][i6]  =  0;  //Initialize  Alp 


for  (i5=0;  i5<NumVisibleSV;  i5++) 

{ 

Alp[i5][0]  =  Dx[i5]; 
Alp[i5][l]  =  Dy[i5]; 
Alp[i5][2]  =  Dz[i5]; 
Alp[i5][3]  =  Dt[i5]; 


//Transpose  Alp  to  get  Brv 
int  i7,  i8; 

double  Brv[4][31]; 


for  (i7=0;  i7<4;  i7++) 

{ 


for  (i8=0;  i8<31;  i8++) 

{ 


Alp[i5][i6]  =  0;  //Initialize  Brv 
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for  (i7=0;  i7<4;  i7++) 

{ 

i8  =0; 

while  (i8<NumVisibleSV) 

{ 

Brv[i7][i8]  =  Alp[i8][i7]; 

i8++; 


//Matrix  multiplication  of  Brv  and  Alp 
int  i9,  ilO,  ill; 
double  Chl[4][4]; 

for  (i9=0;  i9<4;  i9++) 

{ 

for  (il(M);ilO<4;ilO++) 

{ 

Chl[i9][il0]  =  0;  //Initialize  Chi 


for  (i9=0;  i9<4;  i9++) 

{ 


for  (i  10=0;  il0<4;  il0++) 

{ 


for  (i  1 1=0;  ill<31;  ill++) 

Chl[i9][il0]  =  Chl[i9] [i  1 0]  +  Brv[i9][il  l]*Alp[il  l][il0]; 


//Inverse  Chi 
Matrix  m(4,4); 
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int  i  12,  i  13; 


for  (i  12=0;  il2<4;  il2++) 

{ 

for  ( i  1 3  =0 ;  il3<4;  il3++) 

{ 

m(il2+l,  il3+l)  =  Chl[il2][il3];  //Assign  Chi  to  matrix,  m 


Dlt 


Matrix  minv(4,4); 

double  det=Matrix:  :inv(m,minv); 

double  Dlt[4]  [4]; 

for  (i  12=0;  il2<4;  il2++) 

{ 

for  (i  13=0;  il3<4;  il3++) 

{ 

Dlt[il2][il3]  =  minv(il2+l,  il3+l);  //Assign  inversed  matrix,  minv,  to 


//Calculate  DOP 

xdop2  =  sqrt(Dlt[0][0]); 

ydop2  =  sqrt(Dlt[  1  ]  [  1  ]); 

zdop2  =  sqrt(Dlt[2][2]); 

tdop2  =  sqrt(Dlt[3][3]); 

hdop2  =  sqrt(Dlt[0][0]  +  Dlt[l][l]); 

pdop2  =  sqrt(Dlt[0][0]  +  Dlt[l][l]  +  Dlt[2][2]); 

gdop2  =  sqrt(Dlt[0][0]  +  Dlt[l][l]  +  Dlt[2][2]  +  Dlt[3][3]); 

//Assign  x  and  y  values  for  graph  plotting 
XDOP_PlotPoints[il4]=hrl;  //assign  x-values 
XDOP_PlotPoints[il4+l]=  xdop2;  //assign  y-values 
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YDOP_PlotPoints[i  14]=hrl ; 
YDOP_PlotPoints[il4+l]=  ydop2; 

VDOP_PlotPoints[i  14]=hrl ; 
VDOP_PlotPoints[il4+l]=  zdop2; 

TDOP_PlotPoints[il4]=hrl; 
TDOP_PlotPoints[il4+l]=  tdop2; 

HDOP_PlotPoints[i  14]=hrl ; 
HDOP_PlotPoints[il4+l]=  hdop2; 

PDOP_PlotPoints[i  1 4]=hrl ; 
PDOP_PlotPoints[il4+l]=  pdop2; 

GDOPPlotPoint  s  [i  14]  =hr  1 ; 
GDOP_PlotPoints[il4+l]=  gdop2; 

hrl++; 

} 

got_data  =  true; 

II End  of  DOP  calculations  for  24hr  period 

// - 

//Plot  Graph  for  HDOP  as  default 

// - 

m_pPlotItems  =  new  double[num_items*2]; 

//Initialize  the  plot  points 
mlpfs  =  NULL; 
m_pPlotItems  =  NULL; 

m_pPlotItems  =  HDOP  PlotPoints;  //Assign  plot  points 
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//Set  the  graph  position  (left-right,  up-down)  and  size  (width,  height) 
m_pGraph->SetGraphSizePos(230,230,0,0); 


m  lpfs  =new  G_FUNCTIONSTRUCT; 
memset(m_lpfs,0,sizeof(G_FUNCTIONSTRUCT)); 
m_lpfs->FuncType=G_PLOTXY ; 

m_lpfs->szGraphTitle=" Variation  in  F1DOP  over  24hr  Period"; 

m_lp  fs->xMax=2  4 ; 

m_lpfs->xMin=0; 

m_lpfs->yMax=5 ; 

m_lpfs->yVImA); 

m_lpfs->ChartType=G_LINECHART; 

m_lpfs->pPlotXYItems=m_pPlotItems; 

m_lpfs->num_PlotXYItems=num_items; 

m_pGraph->DoFunction(m_lpfs); 

//End  of  graph  plotting 

//Set  display  to  indicate  F1DOP  is  checked 

((CButton*)GetDlgItem(IDC_XDOP_RADIO))->SetCheck(0); 

((CButton*)GetDlgItem(IDC_YDOP_RADIO))->SetCheck(0); 

((CButton*)GetDlgItem(IDC_VDOP_RADIO))->SetCheck(0); 

((CButton*)GetDlgItem(IDC_TDOP_RADIO))->SetCheck(0); 

((CButton*)GetDlgItem(IDC_PDOP_RADIO))->SetCheck(0); 

((CButton*)GetDlgItem(IDC_GDOP_RADIO))->SetCheck(0); 

((CButton*)GetDlgItem(IDC_HDOP_RADIO))->SetCheck(  1 ); 


void  CDOP  CalculatorDlg:  :OnBnClickedXdopRadio() 

{ 

if  (got_data  ==  true) 

{ 


//Initialize  the  plot  points 
m_lpfs  =  NULL; 
m_pPlotItems  =  NULL; 
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m_pPlotItems  =  XDOPPlotPoints;  //Assign  plot  points 

//Set  the  graph  position  (left-right,  up-down)  and  size  (width,  height) 
m_pGraph->SetGraphSizePos(230,230,0,0); 

mlpfs  -new  G_FUNCTIONSTRUCT; 
memset(m_lpfs,0,sizeof(G_FUNCTIONSTRUCT)); 
ml  pfs->Fu  ncT  ype=G_PLOTXY ; 

ml  pfs->szGraphTitlc^'' Variation  in  XDOP  over  24hr  Period"; 

m_lp  fs->xMax=2  4 ; 

m_lpfs->xMin=0; 

m_lpfs->yMax=5 ; 

m_lpfs->yMin=0; 

m_lpfs->ChartType=G_LINECHART; 

m_lpfs->pPlotXYItems=m_pPlotItems; 

m_lpfs->num_PlotXYItems=num_items; 

m_pGraph->DoFunction(m_lpfs); 

//End  of  graph  plotting 


else 

{ 

AfxMessageBox((CString)"Error:  No  data  to  plot  graph.  Press  'Compute' 
Button.",  MB  JCONERROR); 


void  CDOP_CalculatorDlg::OnBnClickedYdopRadio() 

{ 

if  (got_data  ==  true) 

{ 

//Initialize  the  plot  points 
mlpfs  =  NULL; 
m_pPlotItems  =  NULL; 

m_pPlotItems  =  YDOP  PlotPoints;  //Assign  plot  points 
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//Set  the  graph  position  (left-right,  up-down)  and  size  (width,  height) 
m_pGraph->SetGraphSizePos(230,230,0,0); 

mlpfs  -new  G_FUNCTIONSTRUCT; 
memset(m_lpfs,0,sizeof(G_FUNCTIONSTRUCT)); 
ml  p  fs->  Fu  n  cT  ype=G_P  LOTX  Y ; 

m_lpfs->szGraphTitle=" Variation  in  YDOP  over  24hr  Period"; 

mlp  fs->xMax=2  4 ; 

m_lpfs->xMinU); 

m_lpfs->yMax=5 ; 

m_lpfs->yMin=0; 

m_lpfs->ChartType=G_LINECHART; 

m_lpfs->pPlotXYItems:=m_pPlotItems; 

m_lpfs->num_PlotXYItems=num_items; 

m_pGraph->DoFunction(m_lpfs); 

//End  of  graph  plotting 


else 

{ 

AfxMessageBox((CString)"Error:  No  data  to  plot  graph.  Press  'Compute' 
Button.",  MB  JCONERROR); 

} 

} 

void  CDOP  CalculatorDlg:  :OnBnClickedVdopRadio() 

{ 

if  (got_data  ==  true) 

{ 

//Initialize  the  plot  points 
mlpfs  =  NULL; 
m_pPlotItems  =  NULL; 

m_pPlotItems  =  VDOP  PlotPoints;  //Assign  plot  points 
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//Set  the  graph  position  (left-right,  up-down)  and  size  (width,  height) 
m_pGraph->SetGraphSizePos(230,230,0,0); 

m  lpfs  -new  G_FUNCTIONSTRUCT; 
memset(m_lpfs,0,sizeof(G_FUNCTIONSTRUCT)); 
m_lpfs->FuncType=G_PLOTXY ; 

m_lpfs->szGraphTitle=" Variation  in  VDOP  over  24hr  Period"; 

m_lp  fs->xMax=2  4 ; 

m_lpfs->xMin=0; 

m_lpfs->yMax=5 ; 

m_lpfs->yMin=0; 

m_lpfs->ChartType=G_LINECHART; 

m_lpfs->pPlotXYItems=m_pPlotItems; 

m_lpfs->num_PlotXYItems=num_items; 

m_pGraph->DoFunction(m_lpfs); 

//End  of  graph  plotting 


else 

{ 

AfxMessageBox((CString)"Error:  No  data  to  plot  graph.  Press  'Compute' 
Button.", MB  JCONERROR); 


void  CDOP  CalculatorDlg:  :OnBnClickedTdopRadio() 

{ 

if  (got_data  ==  true) 

{ 

//Initialize  the  plot  points 
m_lpfs  =  NULL; 
m_pPlotItems  =  NULL; 

m_pPlotItems  =  TDOP  PlotPoints;  //Assign  plot  points 

//Set  the  graph  position  (left-right,  up-down)  and  size  (width,  height) 
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m_pGraph->SetGraphSizePos(230,230,0,0); 


m  lpfs  -new  G  FUN CTION STRUCT ; 
memset(m_lpfs,0,sizeof(G_FUNCTIONSTRUCT)); 
ml  p  fs->  Fu  n  cT  ype=G_P  LOTX  Y ; 

m_lpfs->szGraphTitle=" Variation  in  TDOP  over  24hr  Period"; 

rrilp  fs->xMax=2  4 ; 

m_lpfs->xMinA); 

m_lpfs->yMax=5 ; 

mJpfs->yVIiiwO; 

m_lpfs->ChartType=G_LINECHART; 

m_lpfs->pPlotXYItems=m_pPlotItems; 

m_lpfs->num_PlotXYItems=num_items; 

m_pGraph->DoFunction(m_lpfs); 

//End  of  graph  plotting 


else 

{ 

AfxMessageBox((CString)"Error:  No  data  to  plot  graph.  Press  'Compute' 
Button.",  MB  JCONERROR); 


void  CDOP  CalculatorDlg:  :OnBnClickedl-[dopRadio() 

{ 

if  (got_data  ==  true) 

{ 

//Initialize  the  plot  points 
m_lpfs  =  NULL; 
m_pPlotItems  =  NULL; 

m_pPlotItems  =  HDOP  PlotPoints;  //Assign  plot  points 

//Set  the  graph  position  (left-right,  up-down)  and  size  (width,  height) 
m_pGraph->SetGraphSizePos(230,230,0,0); 
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mlpfs  -new  G  FUN CTION STRUCT ; 
memset(m_lpfs,0,sizeof(G_FUNCTIONSTRUCT)); 
ml  pfs->Fu  ncT  ype=G_PLOTXY ; 

m_lpfs->szGraphTitle=" Variation  in  HDOP  over  24hr  Period"; 

m_lp  fs->xMax=2  4 ; 

m_lpfs->xMin=0; 

m_lpfs->yMax=5 ; 

m_lpfs->yMin=0; 

m_lpfs->ChartType=G_LINECHART; 

m_lpfs->pPlotXYItems=m_pPlotItems; 

m_lpfs->num_PlotXYItems=num_items; 

m_pGraph->DoFunction(m_lpfs); 

//End  of  graph  plotting 


else 

{ 

AfxMessageBox((CString)"Error:  No  data  to  plot  graph.  Press  'Compute' 
Button.  ",MB_ICONERROR); 


void  CDOP  CalculatorDlg:  :OnBnClickedPdopRadio() 

{ 

if  (got_data  ==  true) 

{ 

//Initialize  the  plot  points 
mlpfs  =  NULL; 
m_pPlotItems  =  NULL; 

m_pPlotItems  =  PDOP  PlotPoints;  //Assign  plot  points 

//Set  the  graph  position  (left-right,  up-down)  and  size  (width,  height) 
m_pGraph->SetGraphSizePos(230,230,0,0); 
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mlpfs  -new  G  FUN CTION STRUCT ; 
memset(m_lpfs,0,sizeof(G_FUNCTION STRUCT)); 
ml  pfs->  Fu  n  cT  ype=G_P  LOTX  Y ; 

m_lpfs->szGraphTitle=" Variation  in  PDOP  over  24hr  Period"; 

mlp  fs->xMax=2  4 ; 

m_lpfs->xVIinU); 

m_lpfs->yMax=5 ; 

m_lpfs->yMin=0; 

m_lpfs->ChaitType=G_LINECHART; 

m_lpfs->pPlotXYItems=m_pPlotItems; 

m_lpfs->num_PlotXYItems=num_items; 

m_pGraph->DoFunction(m_lpfs); 

//End  of  graph  plotting 


else 

{ 

AfxMessageBox((CString)"Error:  No  data  to  plot  graph.  Press  'Compute' 
Button.  ",MB_ICONERROR); 


void  CDOP  CalculatorDlg:  :OnBnClickedGdopRadio() 

{ 

if  (got_data  ==  true) 

{ 

//Initialize  the  plot  points 
mlpfs  =  NULL; 
m_pPlotItems  =  NULL; 

m_pPlotItems  =  GDOP  PlotPoints;  //Assign  plot  points 

//Set  the  graph  position  (left-right,  up-down)  and  size  (width,  height) 
m_pGraph->SetGraphSizePos(230,230,0,0); 

m  lpfs  =new  G_FUNCTIONSTRUCT; 
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memset(m_lpfs,0,sizeof(G_FUNCTIONSTRUCT)); 
ml  pfs->  Fu  n  cT  ype=G_P  LOTX  Y ; 

m_lpfs->szGraphTitle=" Variation  in  GDOP  over  24hr  Period"; 

m_lp  fs->xMax=2  4 ; 

m_lpfs->xMin=0; 

m_lpfs->yMax=5 ; 

m_lpfs->ylVlin=0; 

m_lpfs->ChartType=G_LINECHART; 

m_lpfs->pPlotXYItems=m_pPlotItems; 

m_lpfs->num_PlotXYItems=num_items; 

m_pGraph->DoFunction(m_lpfs); 

//End  of  graph  plotting 


else 

{ 

AfxMessageBox((CString)"Error:  No  data  to  plot  graph.  Press  'Compute' 
Button.",  MB  JCONERROR); 

} 


B.  DOPC  ALCULAT  ORDLG.H 

//  DOPCalculatorDlg.h  :  header  file 

#include  "Graph. h" 

#pragma  once 

//  CDOPCalculatorDlg  dialog 

class  CDOP  CalculatorDlg  :  public  CDialog 

{ 

//  Construction 
public: 

double  *m_pPlotItems; 

LPG  FUN CTION STRUCT  m  lpfs; 

CGraph  *m_pGraph; 

CDOP_CalculatorDlg(CWnd*  pParent  =  NULL);  //  standard  constructor 
//  Dialog  Data 

enum  { IDD  =  IDDDOPCALCULATORDIALOG  }; 
protected: 

virtual  void  DoDataExchange(CDataExchange*  pDX);  //  DDX/DDV  support 

//  Implementation 
protected: 

HICON  mhlcon; 

//  Generated  message  map  functions 
virtual  BOOL  OnInitDialog(); 

afx  msg  void  OnSysCommand(UINT  nID,  LPARAM  IParam); 

afx  msg  void  OnPaint(); 

afx  msg  HCURSOR  OnQueryDraglconQ; 

DECLARE  _MESSAGE_MAP() 

public: 

double  yr; 
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double  mth; 
double  day; 
double  hr; 
double  mi; 
double  timezone; 
double  longi; 
double  lat; 
double  alt; 
double  obs; 

BOOL  daylightsaving; 

CString  xdopO; 

C  String  ydopO; 

CString  vdopO; 

CString  tdopO; 

CString  hdopO; 

CString  pdopO; 

CString  gdopO; 

afx  msg  void  OnBnClickedAlmanac(); 
afx  msg  void  OnBnClickedCompute(); 
afx  msg  void  OnBnClickedXdopRadioQ; 
afx  msg  void  OnBnClickedYdopRadioQ; 
afx  msg  void  OnBnClickedVdopRadio(); 
afx  msg  void  OnBnClickedTdopRadioQ; 
afx  msg  void  OnBnClickedHdopRadioQ; 
afx  msg  void  OnBnClickedPdopRadio(); 
afx  msg  void  OnBnClickedGdopRadioQ; 

//Variables  not  automatically  generated  by  Visual  C++ 

CString  almanac  location; 

int  num_items; 

double  daylightsav; 

bool  got  data; 

double  *XDOP_PlotPoints; 
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double  *YDOP_PlotPoints; 
double  *VDOP_PlotPoints; 
double  *TDOP_PlotPoints; 
double  *HDOP_PlotPoints; 
double  *PDOP_PlotPoints; 
double  *GDOP_PlotPoints; 


double  xdop; 
double  ydop; 
double  vdop; 
double  tdop; 
double  hdop; 
double  pdop; 
double  gdop; 
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C.  GRAPH.CPP  [AFTER  20] 

//  Graph. cpp:  implementation  of  the  CGraph  class. 


#include  "stdafx.h" 

#include  <math.h> 

#include  "DOP  Calculator.h" 

#include  "Graph.h" 

#ifdef  DEBUG 

#undef  THISFILE 

static  char  THIS_FILE[]=_FILE_; 

#define  new  DEBUG  NEW 

#endif 

#pragma  warning  (disable:4244  4018  4701) 

CGraph : :  CGraph( ) 

{ 

GraphSetAllDefaults(); 

} 


/////////////////////  Constructor  ////////////////////////////////////// 

/* 

*/ 

//////////////////////////////////////////////////////////////////////// 

CGraph ::CGraph(CWnd  *pParentWnd,  int  xPos,  int  yPos,  int  Width,  int  Fleight,  UINT  colorscheme) 

{ 

//First  set  everything  to  their  default  values 
GraphSetAllDefaultsQ; 

//setup  default  values 
m_p  W  nd=pParentW  nd; 

SetColorScheme(colorscheme); 

//set  graph  position 
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mi  GraphX=xPo  s<0  ?0 :  xPo  s ; 
miGraph  Y  =yPo  s<0?0:yPos; 

//set  graphsize  (0  or  less  means  default) 

miGraphWidth  =  Width  <  G_MINGRAPHWIDTH?G_MINGRAPHWIDTH:Width; 
m_iGraphHeight=Height  <  G_MINGRAPHHEIGHT?G_MINGRAPHHEIGHT:Height; 


} 

////////////////////  Default  Destructor  /////////////////////////// 
/* 


*/ 


/////////////////////////////////////////////////////////////////// 

CGraph::~CGraph() 

{ 


///////////////////////  GraphSetAllDefaults  /////////////////////////////// 

/*  This  function  sets  all  the  member  variables  to  their  default  values  creates  the  default  font,  etc.... 

Because  there  are  so  may  member  variables,  attempting  to  make  constructors  to  cover  all 
possibilites  will  be  a  pain  -  so  all  constructors  should  call  this  function  first,  then  overwrite  the 
default  values  as  required. 

*/ 

/////////////////////////////////////////////////////////////////////////// 

void  CGraph:  :GraphSetAllDefaults() 

{ 

//setup  default  values 
m_pWnd=NULL; 

//create  default  font 

CreateGraphFont((const  wchar_t  *)"Courier",8); 

SetDefaultColorScheme(); 

//m_bAutofit=TRUE; 

m_bShowGrid=TRUE; 

m_bShowTicks=TRUE; 

//set  default  graph  position 
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m_iGraphX=0; 

m_iGraphY=0; 

//set  default  graph  size 

SetGraphSizePos(0,  0,G_MINGRAPHWIDTH,G_MINGRAPHHEIGHT); 
//set  the  axis  scaling 

SetXAxisScale(G_DEFAULTXMIN,G_DEFAULTXMAX); 

SefYAxisScale(G_DEFAULTYMIN,G_DEFAULTYMAX); 

//set  the  legend  and  Title  Texts 
SetDefaultGraphTitle(); 

SetDefaultXLegendQ; 


//other  stuff 

m_pFunctionParams=NULL; 

} 

///////////////////  SetColorScheme  //////////////////////////////////////////// 
/* 


SetColorScheme  set  the  color  scheme  for  the  graph 

Default  colour  scheme  is  the  same  as  as  the  WFUTE  colorscheme  (which  is  really  grey) 


*/ 

/////////////////////////////////////////////////////////////////////////////// 

void  CGraph::SetColorScheme(int  Scheme,  BOOL  bRedraw) 

{ 

/* 

This  sets  up  the  colors  for  various  graph  elements 

*/ 

switch  (Scheme) 

{ 

case  G_DEFAULTSCHEME: 
case  G_WHITESCHEME: 

m_erYTickColor=RGB(0,0,0); 

m_crXTickColor=RGB(0,0,0); 

m_crYLegendTextColoi=RGB(0,0,0); 

m_crXLegendTextColor=RGB(0,0,0); 
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mcrGraphT  itleColor=RGB  ( 0,0,0 ) ; 
m_crGraphPenColor=RGB(0,0,0); 
mcrGraphBkC  olor=RGB  ( 1 92, 1 92, 1 92); 
m_crGridColor=RGB(220,220,220); 
break; 

default: 

break; 


} 

if(bRedraw) 

{ 


PaintGraphQ; 


/////////////////////////////  CreateGraphFont ///////////////////////////////////// 
/* 


CreateGraphFont  -  Creates  a  font  using  the  specified  facename  and  point  size 


*/ 

/////////////////////////////////////////////////////////////////////////////////// 

void  CGraph::CreateGraphFont(CString  FaceName,  UINT  size) 

{ 

//at  this  point  we  may  not  have  a  handle  to  a  window  -  so  to  fill 
//things  like  char-width  &  height  where  we  need  a  dc  we  will  get  a 
//whole  screen  dc 

CFont  *poldfont; 

HDC  hDC; 

CDC  *dc; 

TEXTMETRIC  textmetrics; 

if  (!m_pWnd)//if  we  don't  have  a  holding  window 

{ 

hDC=GetDC(0);//get  a  whole  screen  dc 
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dc=new  CDC; 
dc->Attach(hDC); 

} 

else 

{ 

dc=m_pWnd->GetDC(); 

} 


//detach  old  font  if  any 
mGraphFont.DetachQ; 

//create  the  new  one 

m_GraphFont.CreatePointFont(size*10,FaceName,dc); 

poldfont=dc->SelectObject(&m_GraphFont); 

dc->GetTextMetrics(&textmetrics); 
miGharHe  ight =t  ext  m  e  t  r  i  c  s .  t  m  Height; 
m_iCharWidth=textmetrics.tmAveCharWidth; 

dc->SelectObject(poldfont); 

if(!m_pWnd) 

{ 

dc->Detach(); 

ReleaseDC(0,hDC); 
delete  dc; 

} 

else 

{ 

m_pWnd->ReleaseDC(dc); 

} 


m_iFontSize=size; 
mszF  ontFace^FaceName; 

//we  need  to  rescale  the  graph 
SetXAxisScale(m_dXAxisMin,  m  dXAxisMax); 
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SetYAxisScale(m_dYAxisMin,  mdYAxisMax); 


} 

//////////////////////////  SetGraphSizePos  ///////////////////////////////////////////// 
/* 


Sets  the  size  of  the  entire  graph  size  and  position  with  the  holding  parent  window. 

A  value  of  -1  for  any  parameter  means  don't  change  that  parameter. 

A  value  of  0  for  the  Width  or  Height(or  any  value  between  0  the  minimum  graph  width  or  height) 
means  use  the  default  graph  width  or  height 

*/ 

/////////////////////////////////////////////////////////////////////////////////// 

void  CGraph::SetGraphSizePos(int  xPos,  int  yPos,  int  Width,  int  Height) 

{ 


CRect  rect; 

rcct.  I  ctt — m  i  Graph  X ; 

rcct.top_m_i  Graph  Y ; 

rcct.  right^rect.  Iett-m_i  Graph  Width; 

rect.bottom=rect.top+m_iGraphHeight; 

/* 

If  xPos  or  yPos  <=0  then  position  will  not  be  changed 

*/ 

m  i  Graph  X-xPos  <  0?m_iGraphX:xPos; 
m_iGraphY=yPos  <  0?m_iGraphY:yPos; 

/* 

A  negative  number  or  zero  means  no  change  of  current  width  or  height 
if  a  +ive  size  given  thats  smaller  than  the  default  min  size  for 
height  or  width  then  use  the  default  min  size  for  those  parameters; 

*/ 

Widths  (Width>0  &&  Width  <G_MINGRAPHWIDTH)?G_MINGRAPHWIDTH:Width; 
miGraphWidth  =  Width  <=  0?  miGraphWidth  :Width; 

Height  =  (Height>0  &&  Height  <G_MINGRAPHHEIGHT)?G_MINGRAPHHEIGHT:Height; 
m  iGraphHeight  =  Height  <=  0?  miGraphHeight :  Height; 
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//remember  to  rescale  the  graph 
S  etXAxisS  cale(m_dXAxisMin,m_dXAxisMax) ; 

S  et  Y  AxisS  cale(m_d  Y  AxisMin,m_d  Y  AxisMax) ; 

//clear  old  graph  from  the  screen 
if(m_pWnd) 

{ 

m_pWnd->InvalidateRect(&rect,TRUE); 
m_pWnd->SendMessage(WM_P  AIN'T, 0,0); 

} 


//////////////////////////  CalcTopMargin  /////////////////////////////////// 

/* 

Returns  the  spacing  between  the  top  edge  of  the  graph  and  the  top  of  the  plotting  area. 

*/ 

//////////////////////////////////////////////////////////////////////////// 

UINT  CGraph:  :CalcTopMargin() 

{ 

return  4*m_iCharHeight; 

} 


/////////////////////  CalcBottommargin  ////////////////////////////////// 
/* 


Returns  the  space  between  the  bottom  of  the  graph  and  the  bottom  of  the  actual  plotting  area 

*/ 

////////////////////////////////////////////////////////////////////////// 

UINT  CGraph:  :CalcBottomMargin() 

{ 

return  3*m_iCharHeight; 

} 


///////////////////////  CalcLeftmargin  //////////////////////////////////// 
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Returns  the  margin  between  the  left  of  the  graph  and  the  left  of  the  plotting  area 


/* 


*/ 

//////////////////////////////////////////////////////////////////////////// 


UINT  CGraph:  :CalcLeftMargin() 

{ 

return  4*m_iCharWidth; 

} 


/////////////////////////  CalcRightMargin  //////////////////////////////////// 

/* 

Returns  the  space  between  the  right  side  of  the  plotting  area  and  the  right  side  of  the  graph 

*/ 

////////////////////////////////////////////////////////////////////////////// 

UINT  CGraph:  :CalcRightMargin() 

{ 

return  m  iCharWidth; 

} 


///////////////////////  SetXAxis Scale  ///////////////////////////////////////// 

/* 

This  sets  the  min  and  max  values  of  the  x-axis  (of  the  plotting  area) 

It  also  sets  what  proportion  of  the  x-range  is  denoted  by  one  pixel  a.ka.  PixelsPerX 

It  also  caluclates  where  the  x-origin  (the  x=0)  point  lies  along  the  x-axis  -  if  x  does  not  pass 
through  zero  the  x-origin  gidline  is  shown  at  the  top  or  bottom  of  the  graph. 

(the  x-origin  point  on  screen  is  stored  in  the  m  iOriginX  parameter  as  an  offset  from  the  LHS  of 
the  (total)  graph 

26/2/2005  -  The  autofit  parameter  is  always  TRUE  as  scrolling  is  not  yet  implemented. 

*/ 

////////////////////////////////////////////////////////////////////////////// 

void  CGraph: :SetXAxisScale(double  min,  double  max) 

{ 

//swap  min  and  max  if  they  are  the  wrong  way  around 
double  temp, scale; 
if  (max  <  min) 

{ 

temp  =min; 
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min=max; 

max=temp; 

} 

//if  min  and  max  are  the  same  (especially  if  they  are  both  zero 
//it  can  be  a  problem  -  so  give  them  a  bit  of  room 
if(min==max) 

{ 

max+=  0.1; 

} 


//set  the  member  variables 

m_dXAxisMax=max; 

m_dXAxisMin=min; 

UINT  lmargin=CalcLeftMargin(); 

UINT  rmargin=CalcRightMargin(); 

if  (m_bAutofit)  //Autofit  is  always  true  @  26  Feb  2005 

{ 

temp=max-min;//the  spread  of  the  x-axiz 
scale=(m_iGraphWidth-lmargin-rmargin)/temp;//calc  pixels/x 
m_dPixelsPerX=scale; 

} 

else 

{ 

m_dPixelsPerX=l  ;//l :  1 

} 

//where  would  the  x-origin  be  located? 
if(  (min  <  0)  &&  (max  >0)) 

{ 

m_iOriginX=abs(min)  *  m  dPixelsPerX  +  lmargin; 

} 

else  if  ( (min  <0)  &&  (max  <=0)) 

{ 
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m_iOriginX=(m_iGraphWidth)-rmargin; 

} 

else  if  (min  >=0  &&  (max>=0)) 

{ 

m_iOriginX=lmargin; 


///////////////////////  Set YAxis Scale  /////////////////////////////////////// 

/* 

This  sets  the  min  and  max  values  of  the  y-axis  (of  the  plotting  area) 

It  also  sets  what  proportion  of  the  grapgh  axis  pixels  represents  1 Y 
(PixelsPerY  =plotheightinpixels/yrange) 

It  also  caluclates  where  the  y-origin  (the  y=0)  point  lies  along  the  y-axis  -  if  y  does  not  pass 
through  zero  the  y-origin  gidline  is  shown  at  the  left  or  right  of  the  graph  (this  can  be  overridden 
by  using  the  SetYLineAtLeft()  function. 

(the  y-origin  point  on  screen  is  stored  in  the  m_iOriginY  parameter  as  an  offset  from  the  bottom  of 
the  (total)  graph 

26/2/2005  -  The  autofit  parameter  is  always  TRUE  as  scrolling  is  not  yet  implemented. 

*/ 

///////////////////////////////////////////////////////////////////////////// 

void  CGraph::SetYAxisScale(double  min,  double  max) 

{ 

double  temp, scale; 

//swap  min  and  max  if  they  are  the  wrong  way  around 
if  (max  <  min) 

{ 

temp=min; 

min=max; 

max=temp; 

} 

//if  min  and  max  are  the  same  (especially  if  they  are  both  zero 
//it  can  be  a  problem  -  so  give  them  a  bit  of  room 
if(min==max) 

{ 

max+=  0.1; 
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} 

//set  the  member  variables 
m_dYAxisMin=min; 
m_dY  AxisMax=max; 

//calculate  scaling 

UINT  bmargin=CalcBottomMargin(); 

UINT  tmargin=CalcTopMargin(); 

CRect  dataarea=CalcDataArea(); 

if  (m_bAutofit)//Always  TRUE 

{ 

temp=max-min;//the  spread  of  the  y-axiz 

scale=(m_iGraphHeight-(bmargin+tmargin))/temp;//calc  pixels/x 
m_dPixelsPerY=scale ; 

} 

else 

{ 

m_dPixelsPerY=l  ;//l :  1 

} 


//where  should  the  Y  origin  be? 

if  (min<0  &&  max  >0)//if  Y  passes  through  zero 

{ 

//from  the  bottom  of  the  graph 
m_iOriginY=abs(min)*m_dPixelsPerY+bmargin; 

} 

else  if(min<0  &&  max<=0)//if  Y  values  are  all  negative 

{ 

m_iOriginY=(dataarea.bottom-dataarea.top)+bmargin; 

} 

else  if  (min  >=0  &&  max  >=0)//if  Y  values  are  all  positive 

{ 

m_iOriginY=bmargin; 
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////////////////////////  PaintGraph  ////////////////////////////////////// 
/* 


This  paints  the  entire  graph  on  to  the  holding  window's  client  area; 

It  does  it  in  steps  starting  from  the  background  and  working  forward. 

As  the  graph  is  NOT  a  window  object  in  it's  own  right,  it  uses  the  display  context  of  the  holding 
window.  If  it  has  not  been  given  a  pointer  to  the  holding  window,  it  will  not  paint. 

Any  CGraph  routine  that  paints  to  the  screen,  checks  the  window  pointer  first. 

The  last  thing  to  be  painted  is  the  plotting  of  the  function  data  (if  any) 

*/ 

/////////////////////////////////////////////////////////////////////////// 

void  CGraph:  :PaintGraph() 

{ 

//here  we  draw  the  graph 

//step  1:  Draw  the  surrounding  rectangle 

//for  the  whole  graph 

if  (m_pWnd==NULL) 

{ 

return; 

} 

CRect  rect; 

CPen  pen,  *oldpen; 

CDC  *dc=m_p W nd->GetDC ( ) ; 

//some  useful  calculations 
UINT  lmargin=CalcLeftMargin(); 

U1NT  rmargin=CalcRightMargin(); 

UINT  bmargin=CalcBottomMargin(); 

UINT  tmargin=CalcTopMargin(); 

UINT  G raphbottom=m_ i Graph  Y +m_iGraphHeight; 

UINT  G  ra phrigh  tuii  iG ra p h  X +m_i  G rap h  W  i d tli ; 

//step  2:  color  the  background 
CBmsh  bmsh,*poldbrush; 
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brush.  CreateSolidBmsh(mcrGraphBkColor); 

pen.CreatePen(PS_SOLID,l,m_crGraphBkColor); 

rect.left=m_iGraphX; 

rect. right  =rect.left+m_iGraphWidth; 

rcct.top_m_i  Graph  Y ; 

rect.bottom=rect.top+m_iGraphHeight; 

oldpen=dc->SelectObject(&pen); 

poldbmsh=dc->SelectObject(&brush); 

dc->Rectangle(&rect); 

dc->SelectObject(oldpen); 

dc->SelectObject(poldbmsh); 

pen.Detach(); 

//step  3:  Draw  Grid  if  required 
DrawGrid(); 

//step  4:  Draw  Axes 
//draw  x-axis 

dc->MoveTo(m_iGraphX+lmargin,Graphbottom-m_iOriginY); 

pen.CreatePen(PS_SOLID,l,m_crXTickColor); 

oldpen=dc->SelectObject(&pen); 

dc->LineTo(Graphright-rmargin,Graphbottom-m_iOriginY); 

dc->SelectObject(oldpen); 

pen.Detach(); 

//draw  the  Y  Axis 

pen.CreatePen(PS_SOLID,l,m_crXTickColor); 
oldpen=dc->SelectObject(&pen); 
if  (!m_bYLineAtLeft) 

{ 

//draw  the  Y  Line  so  that  it  intercepts 
//the  x-line  like  crosshairs 

dc->MoveTo(m_iGraphX+m_iOriginX,m_iGraphY+tmargin); 

dc->LineTo(m_iGraphX+m_iOriginX,Graphbottom-bmargin); 


else 
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{ 

//draw  the  Y  Line  at  the  LHS 

dc->MoveTo(m_iGraphX+lmargin,m_iGraphY+tmargin); 

dc->LineTo(m_iGraphX+lmargin,Graphbottom-bmargin); 

} 

dc->SelectObject(oldpen); 

pen.DetachQ; 

//step  5:  draw  ticks 
DrawTicks(); 

//step  6:  Write  Graph  title 
DrawGraphTitle(); 

//step  7:  Write  x-legend 
DrawXLegend(); 

//step  8:  write  the  x  &  y  axes  values 
DrawXAxisNumbers(); 

DrawYAxisNumbers(); 

//Step  9 

//draw  Function 
DrawFunction(); 

//Cleanup 

m_pWnd->ReleaseDC(dc); 

} 

////////////////////  DrawGrid  //////////////////////////////////////////////// 

/* 

The  grid  comprises  two  parts  -  the  rectangle  drawn  around  the  plotting  areaand  the  vertical&horizontal 
gridlines. 

The  outline  rectangle  is  always  drawn,  the  drawing  of  the  gidlines  is  contrlled  by  the  m_bShowGrid 
member  parameter  using  the  ShowGrid  function. 

The  grid  color  is  mcrGridColor. 


////////////////////////////////////////////////////////////////////////////// 

void  CGraph::DrawGrid() 

{ 


if(!m_pWnd) 

{ 

return; 

} 


/* 

Always  draw  the  dataarea  outline  rectangle 

*/ 

CRect  dataarea=CalcDataArea();;//where  the  graph  data  is  actually  drawn 

//need  a  pen  of  colour  m  crGridColor 
CPen  pen,  *poldpen; 

pen.CreatePen(PS_SOLID,l,m_crGridColor); 

CDC  *pdc=m_pWnd->GetDC(); 

//to  make  a  rectangle  outline  we  have  to  use  a  polyline 

//need  an  array  of  points 

//a  recatangle  comprises  5  point 

CPoint  points[5]; 

//topleft 

points[0].x=dataarea.left; 
po  ints  [  0  ]  ,y=data  area .  top ; 

//topright 

po  ints  [  1  ].  x=dataarea .  right ; 
points[l].y=dataarea.top; 

//bottomright 

points[2]  ,x=dataarea.  right; 
points[2]  ,y—dataarca.  bottom; 

//leftbottom 

points[3].x=dataarea.left; 
points]  3  |.y_dataarca. bottom; 
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//back  to  topleft 
points[4].x=dataarea.left; 
points[4]  .y=dataarea.top; 

poldpen=pdc->SelectObject(&pen); 
pdc->Polyline(points,5);//draw  the  outline  rectangle 

/*  Now  check  whether  the  grid  itself  should 
be  shown 

*/ 

if(!m_bShowGrid) 

{ 

pdc->SelectObject(poldpen); 
m_p  Wnd->ReleaseDC(pdc) ; 
return; 

} 


//draw  the  X-axis  gridlines 

//note  x-axis  grid  lines  run  top  tom  bottom 

double  GridSpacing; 

GridSpacing=CalcXAxisGridAndTicks(); 
int  n;  //for  the  loop 

for(n=l ;  n<G_X_NUMTICKSANDGRID;n++) 

{ 

pdc->MoveTo(dataarea.left+GridSpacing*n ,  dataarea.top); 
pdc->LineTo(dataarea.left+GridSpacing*n,dataarea.bottom); 

} 

//do  the  Y  grid  lines 
//note  Y  gridlines  run  left  -  right 
GridSpacing=Calc  YAxisGridAndT  icksQ; 
for(n=l ;  n<G_Y_NUMTICKSANDGRID;n++) 

{ 

pdc->MoveTo(dataareadeft ,  dataarea.top+GridSpacing*n); 
pdc->LineTo(dataarea. right, dataarea.top+GridSpacing*n); 
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pdc->SelectObject(poldpen); 
m_p  W  nd->ReleaseDC  (pdc) ; 


} 

//////////////////////  DrawTicks  /////////////////////////////////////////// 

/* 

DrawTicks  does  two  things  -  it  draws  the  x  &  Y  axis  lines  and  it  also  draws  the  little  'tick'  lines. 

The  axis  lines  are  always  shown  but  the  'ticks'  are  controlled  by  the  m_bShowTicks  member 
(using  the  ShowTicks  function); 

The  length  of  the  ticks  are  set  by  the  #define  in  the  header  file  the  axis  lines  and  ticks  use  the  same 
color  (m  crXTickColor  for  the  x-axis  and  m  crYTickColor  for  the  y-axis) 

*/ 

void  CGraph::DrawTicks() 

{ 

//pretty  much  the  same  as  showing  the  grid  because  the  ticks  and 
//tick  align 
if  (!m_pWnd) 

{ 

return; 

} 

if  (!m_bShowTicks) 

{ 

return; 

} 

CPen  TickPen,  *poldpen; 

CRect  dataarea=CalcDataArea(); 

UIN  T  GraphB  o  ttom=m_i  Graph  Y +m_iGraphHeight ; 

CDC  *pdc=m_pWnd->GetDC(); 

//start  with  the  x-axis 

TickPen.CreatePen(PS_SOLID,l,m_crXTickColor); 

poldpen=pdc->SelectObject(&TickPen); 

double  GridSpacing=CalcXAxisGridAndTicks(); 

//the  ticklines  vertically  straddle  the  x-axis 

//two  problems  though  -  if  the  x-line  is  at  or  very  close  to  the 

122 


//top  or  bottom  of  the  dataarea 

UINT  xtoptick=((GraphBottom-m_iOriginY)-dataarea.top  < 
G_TICKLENGTH/2)?(GraphBottom-m_iOriginY)-dataarea.top:G_TICKLENGTH/2; 

UINT  xbottick=(dataarea.bottom-(GraphBottom-m_iOriginY) 
<G_TICKLENGTH/2)?dataarea.bottom-(GraphBottom-m_iOriginY):G_TICKLENGTH/2; 

int  n; 

for(n=  1  ;n<G_X_NUMTICKS  ANDGRID  ;n++) 

{ 

//loop  and  do  the  ticks 
//topticks 

pdc->MoveTo(dataarea.left+n*GridSpacing,GraphBottom-m_iOriginY); 
pdc->LineTo(dataarea.left+n*GridSpacing,GraphBottom-m_iOriginY-xtoptick); 
//bottom  ticks 

pdc->LineTo(dataarea.left+n*GridSpacing,GraphBottom-m_iOriginY+xbottick); 


//now  do  the  x  axis  ticks 
pdc->SelectObject(poldpen); 

TickPen.Detach(); 

TickPen.CreatePen(PS_SOLID,l,m_crYTickColor); 

poldpen=pdc->SelectObject(&TickPen); 

GridSpacing=CalcYAxisGridAndTicks(); 

//the  tick  horizontally  straddle  the  the  Y  axis 
//some  problems  though  -  if  the  y-axis  is  at  or  very  close  to  the 
//left  or  right  of  the  data  area  or  if  the  Y -line  memeber  is  set  to  left 
//handside 

UINT  ylefttick=(  (m_iGraphX+m_iOriginX)-dataarea.left  <G_TICKLEN GTH/2 )? 
(m_iGraphX+m_iOriginX)-dataarea.left:G_TICKLENGTH/2; 

UINT  yrighttick=(  dataarea.right-(m_iGraphX+m_iOriginX)  <G_TICKLEN GTH/2 )? 
dataarea. right-(m_iGraphX+m_iOriginX):G_TICKLENGTH/2; 

//check  for  the  special  case  where  the  y-axis  has  been  forced  to  the  left 

if(m_bYLineAtLeft) 

{ 

ylefttick=0; 

} 
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int  x,y; 

if(m_bYLineAtLeft) 

{ 

x=dataarea.left; 
y—dataarca.  bottom; 

} 

else 

{ 

x=m_iGraphX+m_iOriginX; 

y=dataarea.bottom; 

} 

for(n=l;n<G_Y_NUMTICKSANDGRID;n++) 

{ 

//Loop  and  do  the  y  axis  ticks 

pdc->MoveTo(x,y-(n*GridSpacing)); 

//do  left  side  tick 

pdc->LineTo(x-ylefttick,y-(n*GridSpacing)); 
//do  rightside  tick 

pdc->LineTo(x+yrighttick,y-(n*GridSpacing)); 

} 


//cleanup 

pdc->SelectObject(poldpen); 
m_p  Wnd->ReleaseDC(pdc) ; 

} 


/////////////  CalcDataArea  //////////////////////////////////////////////////// 

/* 

This  function  calculates  the  actual  plotting  area  of  the  graph  this  is  the  graph  area  minus  the 
top,bottom,left  &  right  margins 

Returns:  CRect  with  the  plotting  area  (in  client  area  co-ords) 

*/ 
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/////////////////////////////////////////////////////////////////////////////// 

CRect  CGraph::CalcDataArea() 

{ 

CRect  dataarea; 

dataarea. left=m_iGraphX+CalcLeftMargin(); 

dataarea. right=m_iGraphX+m_iGraphWidth-CalcRightMargin(); 

dataarea. top=m_iGraphY+CalcTopMargin(); 

dataarea.  bottom—m_i  Graph  Y+m_i  Graph  Heigh  t-CalcBottomMargin( ); 
return  dataarea; 

} 

/////////////////////  CalcXAxisGridAndTicks  //////////////////////////////// 

/* 

This  calculates  the  positions  of  the  vertical  gridlines  of  the  plot  area. 

This  is  also  used  for  x-axis  ticks  as  the  ticks  line  up  with  the  gridlines. 

How  many  they  are  is  determined  by  the  G  X  NUMTICKSANDGRID  define  in  the  header  file 
Return:  a  double  denoting  the  x-axis  grid  spacing 

*/ 

/////////////////////////////////////////////////////////////////////////// 

double  CGraph: :  CalcXAxisGridAndT icksQ 

{ 

//the  placing  of  the  ticks  co-incide  with  gridlines 
CRect  rect=CalcDataArea(); 

return  ((double)rect.right-(double)rect.left)/(double)G_X_NUMTICKSANDGRID; 


///////////////////  Calc Y AxisGridAndT icks  ////////////////////////////////// 

/* 

This  calculates  the  positions  of  the  horizontal  gridlines  of  the  plot  area. 

This  is  also  used  for  y-axis  ticks  as  the  ticks  line  up  with  the  gridlines. 

How  many  they  are  is  determined  by  the  G  Y  NUMTICKSANDGRID  define  in  the  header  file 
Return:  A  double  denoting  the  horizontal  gridline  spacing 

*/ 

////////////////////////////////////////////////////////////////////////////// 

double  CGraph: :CalcYAxisGridAndTicks() 
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{ 

//the  placing  of  the  ticks  co-incide  with  gridlines 
CRect  rect=Ca  I  c  Data  A  rca( ) ; 

return  ((double)rect.bottom-(double)rect.top)/(double)G_Y_NUMTICKSANDGRID; 

} 


/////////////////////  DrawGraphTitle  //////////////////////////////////// 

/* 

This  draws  the  Graph  title  string  in  the  color  m_crGraphTitleColor 
This  title  is  centred  abobe  the  plotting  area. 

*/ 

//////////////////////////////////////////////////////////////////////// 

void  CGraph::DrawGraphTitle() 

{ 

//The  graph  title  is  drawn  one  character  line  down 
//centered  left  right  between  the  left  and  right  margins 
if  (m_pWnd==NULL) 

{ 

return; 

} 

UINT  lmargin=CalcLeftMargin(); 

UINT  rmargin=CalcRightMargin(); 

CRect  rect; 

rect.left=m_iGraphX+lmargin; 

rect.  top~m_i  Graph  Y+m_iCharHeight; 

rect.right=m_iGraphX+m_iGraphWidth-rmargin; 

rect.bottom=rect.top+m_iCharHeight; 

//draw  the  title  using  the  specified  colorscheme 

//using  the  graph  font 

CDC  *pdc=m_pWnd->GetDC(); 

CFont  *poldfont; 

//note  we  must  clear  of  any  old  stuff  crap  fom  this  area 
CBmsh  brush; 


126 


brush.  CreateSolidBrush(mcrGraphBkColor); 

//or  bottom  of  a  given  rect 
pdc->FillRect(&rect,&brash); 

pdc->S  etBkMode(TRAN  SP  ARENT ) ; 

pdc->SetTextColor(m_crGraphTitleColor); 

poldfont=pdc->SelectObject(&m_GraphFont); 

pdc->DrawText(m_szGraphTitle,&rect,DT_CENTER|DT_END_ELLIPSIS); 

//cleanup 

pdc->SelectObject(poldfont); 

m_pWnd->ReleaseDC(pdc); 


/////////////////////  DrawXLegend  ////////////////////////////////////////////// 

/* 

The  Xaxis  legend  is  drawn  below  the  plotting  area  below  the  x-axis  scale  numbers 

*/ 

/////////////////////////////////////////////////////////////////////////////// 

void  CGraph::DrawXLegend() 

{ 

//The  x  legend  is  drawn  below  the  data  area 
if  (m_pWnd==NULL) 

{ 

return; 

} 

UINT  lmargin=CalcLeftMargin(); 

U1NT  rmargin=CalcRightMargin(); 

CRect  rect; 

rect.left=m_iGraphX+lmargin; 

rect.right=m_iGraphX+m_iGraphWidth-rmargin; 

rect.  top_m_iGraph  Y+m_i  Graph  Heigh  t-2*m_iCharHeight- 1 ; 

rect.bottom=rect.top+m_iCharHeight; 

//just  2  b  safe  move  the  rect  down  a  bit  so  we  don't 
//interfere  with  the  yaxis  numbers 
rect.top+=l; 
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rect.bottom+=l; 


//draw  the  title  using  the  specified  colorscheme 

//using  the  graph  font 

CDC  *pdc=m_pWnd->GetDC(); 

CFont  *poldfont; 

//note  we  must  clear  of  any  old  stuff  crap  fom  this  area 
CBrush  brush; 

brush.  CreateSolidBrush(m_crGraphBkColor); 

rect.lnflateRect(0,0,l,l);//bcause  fdlrect  does  not  go  right  to  the  right 

//or  bottom  of  a  given  rect 

pdc->FillRect(&rect,&brush); 

rect. DeflateRect(0, 0,1,1); 

pdc->S  etBkMode(TRAN  SPARENT ) ; 
pdc->SetT  extColor(m_crXLegendT  extColor); 
poldfont=pdc->SelectObject(&m_GraphFont); 

pdc->DrawText(m_szXLegendText,&rect,DT_CENTER|DT_END_ELLIPSIS); 

//cleanup 

pdc->SelectObject(poldfont); 
m_p  W  nd->ReleaseDC  (pdc) ; 


///////////////////// DrawXAxisNumbers  ////////////////////////////////////// 

/* 

The  X  axis  scale  numbers  are  drawm  directly  below  the  plot  area. 
Three  numbers  are  drawn,  min,  middle  and  max  scaling 

*/ 

//////////////////////////////////////////////////////////////////////////// 

void  CGraph::DrawXAxisNumbers() 

{ 

//we  draw  three  sets  of  numbers  xmin,  xmiddle  and  xmax 
//we  will  limit  them  to  7  digits 
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// 

if  (!m_pWnd) 

{ 

return; 

} 

CString  astring; 

CRect  rect,dataarea; 

CFont  *poldfont; 

CDC  *pdc; 

pdc=m_p  W  nd->GetDC  () ; 
poldfont=pdc->SelectObject(&m_GraphFont); 
pdc->SetTextColor(m_crXLegendTextColor); 
pdc->S  etBkMode(TRAN  SP  ARENT ) ; 

dataarea=CalcDataArea(); 

//do  the  left  side  (min)  left  aligned 
rect.left~dataarca.left; 

rect.  top~m_iGraph  Y+m_i  Graph  Height-3*m_iCharHeight; 

reGt.right=rect.left+8*m_iCharWidth; 

rect.bottom=rect.top+m_iCharFIeight; 

rect.top+=l; 

rect.bottom+=l; 

//clear  any  old  text 
CBrush  brush; 

brush.  CreateSolidBrush(mcrGraphBkColor); 
pdc->FillRect(&rect,&brush); 

//format  and  print  the  number 

astring.  Format((CString)"%.4g",m_dXAxisMin); 

pdc->DrawText(astring,&rect,DT_NOCLIP|DT_LEFT); 

//do  the  half  way  point  -  centre  aligned 

rect.left=dataarea.left+((dataarea.right-dataarea.left)/2); 

rect.left=rect.left-4*m_iCharWidth; 

rect.right=rect.left+8*m_iCharWidth; 

rect.top+=l; 

rect.bottom+=l; 


129 


pdc->FillRect(&rect,&brush); 

astring.Format((CString)"%.4g",(m_dXAxisMin+m_dXAxisMax)/2); 

pdc->DrawText(astring,&rect,DT_NOCLIP|DT_CENTER); 

//now  do  the  righthand  side  (max);  right  aligned 

rect.left=dataarea.right-8*m_iCharWidth; 

re  ct .  right=dataarea  .right ; 

rect.top+=l; 

rect.bottom+=l;; 

pdc->FillRect(&rect,&bmsh); 

astring.Format((CString)"%.4g",m_dXAxisMax); 

pdc->DrawText(astring,&rect,DT_NOCLIP|DT_RIGF[T); 

//cleanup 

pdc->SelectObject(poldfont); 

m_pWnd->ReleaseDC(pdc); 


} 

lllllllllllllllllllimxwYAxmxxrnbexsIlllllllllllllllllllllllllllllllllllU 

/* 

The  Y  axis  scale  is  drawn  on  the  LF1S  of  the  plot  area. 

Only  two  numbers  are  drawn  -  min  and  max  to  allow  for  the  Y-axis  legend. 

*/ 

///////////////////////////////////////////////////////////////////////////// 

void  CGraph::DrawYAxisNumbers() 

{ 

//we  will  only  do  two  sets  of  numbers  min  and  max  because 
//if  we  print  the  halfway  point  it  will  cross  the  yaxis  label 
if  (!m_pWnd) 

{ 

return; 

} 

CString  astring; 

CRect  rect.dataarca; 

CFont  *poldfont; 
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CDC  *pdc; 


pdc=m_p  W  nd->GetDC  () ; 
poldfont=pdc->SelectObject(&m_GraphFont); 
pdc->SetTextColor(m_crYLegendTextColor); 
pdc->S  etBkMode(TRAN  SPARENT ) ; 

dataarea=CalcDataArea(); 

//do  the  max  value 

rect.left=m_iGraphX-9*m_iCharWidth; 

rect.top=dataarea.top-m_iCharWidth; 

rcct.righ tercet. I eft+1 2*m_iCharWidth;//allow  for  characters 
rect.bottom=rect.top+m_iCharHeight; 

//clear  any  old  text 
CBrush  brush; 

brush.  CreateSolidBrush(mcrGraphBkColor); 

//format  and  print  the  number 
astring.Format((CString)"%.4g",m_dYAxisMax); 
pdc->DrawT  ext(astring,  &rect,DT_RIGFlT ) ; 

//do  the  4/5  value 

rect.top=dataarea.top+2*m_iCharHeight; 

rect.bottom^rect.top+miCharFIeight; 

astring.Format((CString)"%.4g",(0.8*m_dYAxisMax)); 

pdc->DrawText(astring,&rect,DT_RIGF[T); 

//do  the  3/5  value 

rect.top=dataarea.top+4.5*m_iCharHeight; 

rect.bottom^rect.top+miCharFIeight; 

astring.Format((CString)"%.4g",(0.6*m_dYAxisMax)); 

pdc->DrawText(astring,&rect,DT_RIGl-[T); 

//do  the  2/5  value 

rect.top=dataarea.top+7*m_iCharHeight; 

rect.bottom^rect.top+miCharFleight; 
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astring.Format((CString)"%.4g",(0.4*m_dYAxisMax)); 

pdc->DrawText(astring,&rect,DT_RIGHT); 

//do  the  1/5  value 

rect.top=dataarea.top+9.5*m_iCharHeight; 

rect.bottom=rect.top+m_iCharHeight; 

astring.Format((CString)"%.4g",(0.2*m_dYAxisMax)); 

pdc->DrawText(astring,&rect,DT_RIGFlT); 

//do  the  bottom  -  (min) 

rect.top=dataarea.bottom-0.5*m_iCharHeight; 

rect.bottom=rect.top+m_iCharFIeight; 

astring.Format((CString)"%.4g",(m_dYAxisMin)); 

pdc->DrawText(astring,&rect,DT_RIGF[T); 

//cleanup 

pdc->SelectObject(poldfont); 
m_p  W  nd->ReleaseDC  (pdc) ; 


////////////////// SetdefaultColorScheme  ////////////////////////////////// 
/* 


sets  the  colours  of  various  bits  back  to  the  default  colours 


*/ 

/////////////////////////////////////////////////////////////////////////// 

void  CGraph: :  SetDefaultColorSchemeQ 

{ 

SetColorScheme(G_DEFAULTSCHEME); 

} 


//////////////////////  SetDefaultGraphTitl e  ////////////////////////////// 
/* 

Sets  the  default  graph  title 

*/ 
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Illllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllll 

void  CGraph::SetDefaultGraphTitle() 

{ 

m_szGraphTitle="Variation  in  DOP  over  24hr  Period" 


//////////// SetDefaultXLegend ///////////////////////////////////////////// 
/* 


*/ 


/////////////////////////////////////////////////////////////////////////// 

void  CGraph: :SetDefaultXLegend() 

{ 

m_szXLegendText="Time  in  24hr  Format"; 

} 


/////////////////////////  SetGraphTitle  ///////////////////////////////////// 
/* 

This  sets  the  GraphTitle 
Takes: 

C  String 

*/ 

///////////////////////////////////////////////////////////////////////////// 

void  CGraph: :SetGraphTitle(CString  GraphTitle) 

{ 

mszG ra p h  I  itl  c — G rap hTitle ; 

DrawGraphTitle(); 

} 


//////////////////////// SetXLegendText  /////////////////////////////////// 
/* 

Sets  the  X  axis  legend  text 

*/ 
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Illlllllllllllllllllllllllllllllllllllllllllllllilllllllllllllllllllllllllll 

void  CGraph::SetXLegendText(CString  XText) 

{ 

mszXLegendT  ext=XT  ext; 
DrawXLegend(); 

} 


/////////////////////  SetYLineAtLeft  //////////////////////////////////////// 

/* 

The  Y  axis  line  can  be  forced  to  the  LHS  of  the  plot  area  using  this  function 

*/ 

///////////////////////////////////////////////////////////////////////////// 

void  CGraph::SetYLineAtLeft(BOOL  AtLeft) 

{ 

BOOL  bprevious=m_bYLineAtLeft; 
m_bYLineAtLeft= AtLeft; 

//if  there  is  a  change  in  the  Y  line  position  then  we  will  have  to 

//redraw  the  graph 

if  (m  bYLineAtLeft  !=  bprevious) 

{ 

PaintGraphQ; 


//////////////////////////  ShowGrid  ///////////////////////////////////////// 

/* 

This  switches  the  grid  on  or  off  as  set  by  the  bShow  parameter 
The  graph  is  repainted  to  match 

*/ 

///////////////////////////////////////////////////////////////////////////// 

void  CGraph::ShowGrid(BOOL  bShow) 

{ 

//this  is  public  function 

//show  the  graph  grid  if  bShow==TRUE 
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//or  vice-versa 

BOOL  bprevious=m_bShowGrid; 
m_bShowGrid=bShow; 

//if  there  is  a  change  then  repaint 
if  (m_bShowGrid  !=bprevious) 

{ 

PaintGraph(); 

} 

} 


///////////////////  ShowTicks  //////////////////////////////////////////// 

/* 

This  switches  the  x  and  y  axis  'ticks'  on  or  off  as  set  by  the 
bShow  parameter 

*/ 

/////////////////////////////////////////////////////////////////////////// 

void  CGraph::ShowTicks(BOOL  bShow) 

{ 

//this  is  public  function 

//show  the  graph  grid  if  bShow==TRUE 

//or  vice-versa 

BOOL  bprevious=m_bShowTicks; 
m_bShowTicks=bShow; 

//if  there  is  a  change  then  repaint 
if  (m_bShowTicks!=bprevious) 

{ 

PaintGraphQ; 

} 


IIIIIIIIIIIIIIIIIIIIIIConv^oGravhCoorfelllllllllllllllllllllllllllllllllll 

/* 

ConvertToGraphCoords  (double  x,  double  y)  will  take  the  result  of  some  calculation  as  given  by  x  and  y 
and  return  where  they  should  be  plotted  on  the  graph. 
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As  we  are  can  only  plot  a  whole  pixel  the  return  value  is  of  type  LONG  (as  apposed  to  double). 


The  Y-pixel  is  in  the  HIWORD  and  x-pixel  in  the  LOWORD  of  the  return. 

Note  that  depending  on  the  scale  of  the  graph,  the  return  cords  of  a  single  pixel  could  be  one  of  many. 

For  example:  if  the  x  axis  is  400  pixels  wide,  but  is  scaled  to  represent  1000,  then  each  pixel  represent  2.5 
in  the  real  world. 

So  to  pixel  398  represent  the  real  world  values  of  995  to  996.5  inc. 


*/ 

/////////////////////////////////////////////////////////////////////////////// 

LONG  CGraph::ConvertToGraphCoords(double  x,  double  y) 

{ 


LONG  result  =-l; 

//to  be  plottable  on  the  graph  the  given  x-value  must  be  between 
//x-min  and  x-max 

if(x  <  m_dXAxisMin  ||  x  >  m_dXAxisMax) 

{ 

return  result; 

} 

if(  y  <m_dYAxisMin  ||  y  >  m_dYAxisMax) 

{ 

return  result; 

} 


//calc  the  abs  difference  between  Xmin  and  x; 
double  xdif  =  abs(m_dXAxisMin  -  x); 

//calc  the  abs  difference  between  Ymin  and  y; 
double  ydif  =  abs(m_dYAxisMin  -  y); 

//find  the  dataarea 
CRect  rect=CalcDataArea(); 

int  xpos=rect.left+(xdif*m_dPixelsPerX);  //from  left 
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int  ypos=rect.bottom-(ydif*m_dPixelsPerY);  //from  bottom 
result=MAKELONG(xpos,ypos); 

return  result; 


/////////////////////////////  DoFunction  ////////////////////////////////////// 
/* 


The  user  fills  in  a  GFUNCTIONSTRUCT  relevant  to  a  function  and  passes 
a  pointer  to  it  to  this  function. 

The  pointer  to  this  G  FUNCTIONSTRUCT  is  saved  in  a  member  variable 
This  function  does  some  preliminary  stuff  and  if  there  are  no 
obvious  problems  it  then  calls  the  PaintGraph  function. 


Returns: 

FALSE  if  there  are  no  problems 


*/ 

/////////////////////////////////////////////////////////////////////////////// 

BOOL  CGraph::DoFunction(G_FUNCTIONSTRUCT  *pFunctionParams) 

{ 

//do  some  checks  first  to  see  if  good  data  has  been  passed 

if(  (pFunctionParams->pPlotXYItems==NULL)  &&  (pFunctionParams- 
>FuncType=G_PLOTXY)) 

{ 

return  FALSE; 

} 

m_pF  unctionParams=pFunctionParams ; 

//Set  chart  title  and  other  text  items 
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SetGraphTitle((CString)m_pFunctionParams->szGraphTitle); 


////////////////////AutoScaling  X  axis 
/* 

Precautions: 

For  the  plotdeviationpercent  and  plotdeviationabsolute 

the  axis  min  should  be  0  (having  any  other  values  makes  no  sense) 

*/ 

double  xminimum; 

xminimum=pFunctionParams->xMin; 

SetXAxisScale(xminimum,pFunctionParams->xMax); 

//set  the  Y  axis  scale 

SetYAxisScale(pFunctionParams->yMin,pFunctionParams->yMax); 

PaintGraph(); 
return  TRUE; 


/////////////////////////  ClearFunction  //////////////////////////////// 
/* 


This  resets  the  GFUNCTIONSTRUCT  member  pointer 
The  graph  is  repainteed  (cleared) 

*/ 

//////////////////////////////////////////////////////////////////////// 

void  CGraph::ClearFunction() 

{ 

m_pFunctionParams=NULL;  //reset  the  pointer 

//Clear  the  graph 

PaintGraph(); 

} 
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IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIITimsNYvmchon  ///////////////////////////////// 

/* 

//  This  is  called  from  within  the  PaintGraph  routine  to  draw 
//  the  actual  function  onto  the  graph. 

//  This  functions  just  switches  the  FunctionType  member  specified  in  the 

//  G  FUNCTIONSTRUCT  and  calls  the  appropriate  routine. 

*/ 

/////////////////////////////////////////////////////////////////////////////// 

void  CGraph::DrawFunction() 

{ 

if(!m_pWnd) 

{ 

return; 

} 


if(!m_pFunctionParams) 

{ 

return; 

} 

switch  (m_pFunctionParams->FuncType) 

{ 

case  GPLOTXY : 

DoPlotXY(); 

break; 

default: 

break; 

} 


////////////////////////  DoPlotXY  //////////////////////////////////////////////// 

/* 
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This  has  many  similarities  with  the  other  functions  -  however 
there  is  no  need  to  calculate  Y  in  the  sane  way  as  in  the  other  function 
as  it  is  given. 

a  particular  x  point  on  the  graph  is  also  given  -  we  place  the 
given  y  value  at  the  given  x  point 

Other: 

a.  We  will  not  bother  plotting  if  x  <  the  graph  scale  XMin  or  >  XMax 

b.  Only  dot,  bar  and  line  will  be  acceptable  for  the  chart  type.  Any  other 
will  default  to  bar.  However  be  aware  that  line  is  only  suitable  if  the 
x-values  are  in  sequence. 

c.  The  usual  y  constraints  apply 


*/ 


/////////////////////////////////////////////////////////////////////////////////// 

void  CGraph:  :DoPlotXY () 

{ 

UINT  prevx=0; 

UINT  prevy=0; 

BOOL  firstpoint=TRUE; 

LONG  result; 

UINT  xstart; 

UINT  ystart; 

double  xperpixel=  1/mdPixelsPerX; 
double  yperpixel=l/m_dPixelsPerY ; 

result=ConvertToGraphCoords(m_dXAxisMin,m_dYAxisMin); 
xstart=LOWORD(result);//the  left  hand  side  of  the  graph  plot  area  on  screen 
ystart=HIWORD(result);//shouldbe  the  bottom  of  the  graph  plot  area  on  screen 

for  (UINT  count  =0;  count  <  m_pFunctionParams->num  PlotXYItems*2;  count+=2) 

{ 
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double  x=m_pFunctionParams->pPlotXYItems[count]; 
double  y=m_pFunctionParams->pPlotXYItems[count+l]; 


//if  x  is  off  scale  -  don't  bother 

if  ( (x  <  m_dXAxisMin)  ||  (x  >  m_dXAxisMax)) 

{ 

continue;  //NEXT  !!!!!! 

} 


UINT  pixelx=  xstart+(x-m_dXAxisMin)/xperpixel; 
//y=ConstrainY  (y); 

UINT  pixely=ystart-(y-m_dYAxisMin)/yperpixel; 

if(firstpoint) 

{ 

prevx=pixelx; 

prevy=pixely; 

firstpoint=FALSE; 

} 

PlotPoints(pixelx,pixely,prevx,prevy); 

//current  point  becomes  previous  point 

prevx=pixelx; 

prevy=pixely; 


///////////////// PlotPoints  ///////////////////////////////////////////// 
/* 


Each  Function  e.g.  DoSineX,DoPlotXY,  etc,  calls  this  function  as 
they  calculate  each  point  so  that  each  point  can  be  drawn  on  the 
Takes: 


141 


U1NT  x,  UNIT  y  -  the  graph  co-ord  of  the  point  just  calculeted 
(current  point). 

UINT  prevx,  UINT  prevy  -  the  co-ords  of  the  previous  point 
This  routine  checks  what  type  of  plot  (line,  dot,  or  bar)  is 
required  and  calls  the  appropriate  routine 

*/ 

/////////////////////////////////////////////////////////////////////////// 

void  CGraph:  :PlotPoints(UINT  x,  UINT  y,  UINT  prevx,  UINT  prevy) 

{ 

//here  we  check  the  chart  type  and  plot  the  points  accordingly 
//we  need  to  constarin  the  Y  values  to  keep  them  within  the 
//plot  area; 

switch(m_pFunctionParams->ChartType) 

{ 

case  GLINECHART : 

{ 

DrawConnectLine(prevx, prevy, x,y); 
break; 

} 

} //SWITCH 

return; 

} 


////////////////////////////  Dra wConnectL  ine  ////////////////////////////////// 
/* 


For  the  line  chat  type,  this  routine  draws  a  line  between  previous 
point  (FROM)  and  current  point  (TO) 


*/ 

/////////////////////////////////////////////////////////////////////////////// 

void  CGraph:  :DrawConnectLine(UINT  FromX,  UINT  FromY,  UINT  ToX,  UINT  ToY) 
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//draws  a  connecting  line  between  to  pixels 
//using  the  graphpen  color 

if(!m_pWnd) 

{ 

return; 

} 

CPen  pen,  *poldpen; 

pen.CreatePen(PS_SOLID,l,m_crGraphPenColor); 

CDC  *pdc=m_pWnd->GetDC(); 

poldpen=pdc->SelectObject(&pen); 

pdc->MoveTo(FromX,FromY); 

pdc->LineTo(ToX,ToY); 

pdc->SelectObject(poldpen); 

m_pWnd->ReleaseDC(pdc); 


D 


GRAPH.H  [AFTER  20] 


//  Graph.h:  interface  for  the  CGraph  class. 

// 

////////////////////////////////////////////////////////////////////// 

#if  !defined(AFX_GRAPH_H_70FB8DF3_88AC_40C5_802C_58621127B9E9_INCLUDED_) 
#define  AFXGRAPHH _ 7 OFB  8DF3  8  8 AC_40C5_802C_5  862 1 1 27B9E9 _ IN CLUDED_ 

#if_MSC_VER>  1000 
#pragma  once 

#endif  //  MSC  VER>  1000 


//  Some  defines 


//colorschemes 

#define  G_DEF  AULT  SCHEME  0 
#define  G  WHITESCHEME  1 


//default  graphsize 

#define  GMIN  GRAPH  WIDTH  300 
#define  G  MINGRAPHHEIGHT  200 


//default  axies  scaling 
#define  GDEF AULTXMIN  0 
#define  G  DEF AULTXMAX  24 
#define  GDEF  AULT  YMIN  0 
#define  GDEF  AULT  YMAX  5 

//miscellaneous 

#define  G  X  NUMTICKSANDGRID  12  //how  may  parts  the  dataarea  is  divided 
#define  G  Y  NUMTICKSANDGRID  5  //how  may  parts  the  dataarea  is  divided 
#define  G  TICKLENGTH  10  //size  of  those  little  ticks  on  the  axes 

/* 

Function  related  defines  and  stuff 
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*/ 

//for  the  builtin  functions 
#define  NUMFUNCTIONS  1 
//Plot  function 
#define  G_PLOTXY  1 

//plot  type  bar, line,  etc... 

#define  NUMCHARTTYPES  1 

#define  G  LINECHART  1  //  each  point  is  a  line  drawn  from  each  x-y  point  to  the  next 

//some  structures  for  passing  data 
typedef  struct 

{ 

UINT  FuncType; 

U1NT  ChartType;//line 
double  xMin; 
double  xMax; 
double  yMin; 
double  yMax; 
char  *szGraphTitle; 
char  *szXLegend; 
double  *pPlotXYItems; 
int  numPlotXYItems; 

}G_FUNCTIONSTRUCT,  *LPG_FUNCTIONSTRUCT; 


class  CGraph 

{ 

public: 

void  ShowTicks(BOOL  bShow); 
void  ClearFunction(void); 

BOOL  D oFunction(G_FUN CTION STRU CT  *pFunctionParams); 
void  ShowGrid(BOOL  bShow); 
void  SetYLineAtLeft(BOOL  AtLeft); 
void  GraphSetAllDefaultsQ; 
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void  SetXLegendText(CString  XText); 
void  SetGraphTitle(CString  GraphTitle); 
void  PaintGraph(void); 

CGraph(CWnd  *pParentWnd,int  xPos=0,  int  yPos=0,  int  Width  =0,  int  Height=0,  UINT 
colorscheme=G_DEFAULTSCHEME); 

void  SetYAxisScale(double  min,  double  max); 

void  SetXAxisScale(double  min,double  max); 

void  SetGraphSizePos(int  xPos,  int  yPos,  int  Width,  int  Height); 

void  SetColorScheme(int  Scheme,  BOOL  bRedraw=FALSE); 

void  CreateGraphFont(C String  FaceName,UlNT  size); 

CGraph(); 

virtual  ~CGraph(); 


private: 

LONG  ConvertToGraphCoords(double  x,  double  y); 
void  DoPlotXY(); 

void  DrawC onnectLine(UINT  FromX,  UINT  FromY,  UINT  ToX,  UINT  ToY); 
void  PlotPoints(UINT  x,  UINT  y,  UINT  prevx,  UINT  prevy); 
void  DrawFunction(); 
void  SetDefaultColorScheme(void); 

CString  m_szFunctionNameText; 

COLORREF  m  crFunctionNameColor; 
void  DrawYAxisNumbers(); 
void  DrawXAxisNumbers(void); 
void  DrawTicks(void); 
double  CalcYAxisGridAndTicks(void); 

CRect  CalcDataArea(void); 

double  CalcXAxisGridAndTicks(void); 

void  DrawGrid(void); 

void  DrawXLegend(); 

void  SetDefaultXLegend(void); 

void  SetDefaultGraphTitle(void); 

void  DrawGraphTitleQ; 

UINT  CalcRightMargin(); 

UINT  CalcLeftMargin(); 


146 


U1NT  CalcBottomMarginQ; 

UINT  CalcTopMargin(); 

BOOL  m_bShowTicks;//Ticks  are  the  little  things  on  the  x  &  y  axis 
BOOL  mbShowGrid; 

BOOL  mbAutofit; 

BOOL  mbYLineAtLeft; 

CString  m_szFontFace; 

CString  m  szXLegendText; 

CString  m_szGraphTitle; 

COLORREF  m  crYT ickColor; 

COLORREF  m  crXT ickColor; 

COLORREF  m  crYLegendTextColor; 

COLORREF  m  crXLegendTextColor; 

COLORREF  m  crGraphTitleColor; 

COLORREF  m  crGraphPenColor; 

COLORREF  m  crGraphBkColor; 

COLORREF  m  crGridColor; 
int  m_iFontSize;// 
int  m_iGraphWidth;// 
int  m_iGraphFIeight;// 

int  m_iGraphX;//location  of  the  fraph  within  the  window 

int  m_iGraphY ;//location  of  the  graph  within  the  window 

double  m_dXAxisMin;//  the  start  value  of  X 

double  m_dYAxisMin;//start  value  of  Y 

double  m  dXAxisMax; 

double  m  dYAxisMax; 

CWnd  *m_pWnd;//parent/owner 
//Flelper  calculated  values  -  that  is  to  say  that  these 
//values  are  not  passed  in  to  the  graph 
//they  are  calculated  from  other  given 

CFont  m_GraphFont;//default  font  font  created  from  default  fontface,  &  point  size 

int  m_iCharHeight;//calculated  from  the  font 

int  m_iCharWidth;//calculated  from  the  font 

int  m_iOriginX;//location  of  the  origin  within  the  graph 

int  m_iOriginY ;//location  of  the  origin  within  the  graph 

double  m_dPixelsPerY;//scaling 
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double  m_dPixelsPerX;//scaling 
int  miScrollPosX;// 
int  miScrollPosY;// 

//Data  related  variables 

G  FUNCTIONSTRUCT  *m_pFunctionParams; 

}; 

#endif  // ! defined/ AFX  GRAPH  H _ 70FB8DF3_88AC_40C5_802C_58621 127B9E9_INCLUDEDJ) 
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E.  MATRIX.H  [AFTER  21] 

#include  <assert.h>  //  Defines  the  assert  function. 

class  Matrix  { 
public: 

//  Default  Constructor.  Creates  a  1  by  1  matrix;  sets  value  to  zero. 
Matrix  ()  { 

nRow_  =  1 ;  nCol_  =  1 ; 

data_  =  new  double  [1];  //  Allocate  memory 

set(O.O);  //  Set  value  of  data_[0]  to  0.0 

} 


//  Regular  Constructor.  Creates  an  nR  by  nC  matrix;  sets  values  to  zero. 
//  If  number  of  columns  is  not  specified,  it  is  set  to  1 . 

Matrixfint  nR,  int  nC  =  1)  { 

assert(nR  >  0  &&  nC  >  0);  //  Check  that  nC  and  nR  both  >  0. 
nRow_  =  nR;  nCol_  =  nC; 
data_  =  new  double  [nR*nC];  //  Allocate  memory 
assert(data_  !=  0);  //  Check  that  memory  was  allocated 

set(O.O);  //  Set  values  of  data_[]  to  0.0 

} 


//  Copy  Constructor. 

//  Used  when  a  copy  of  an  object  is  produced 
//  (e.g.,  passing  to  a  function  by  value) 
Matrix(const  Matrix&  mat)  { 
this->copy(mat);  //  Call  private  copy  function. 

} 


//  Destructor.  Called  when  a  Matrix  object  goes  out  of  scope. 
~Matrix()  { 

delete  []  data_;  //  Release  allocated  memory 

} 
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//  Assignment  operator  function. 

//  Overloads  the  equal  sign  operator  to  work  with 

//  Matrix  objects. 

Matrix&  operator=(const  Matrix&  mat)  { 
if(  this  ==  &mat )  return  *this;  //  If  two  sides  equal,  do  nothing, 
delete  []  data_;  //  Delete  data  on  left  hand  side 

this->copy(mat);  //  Copy  right  hand  side  to  l.h.s. 

return  *this; 

} 


//  Simple  "get"  functions.  Return  number  of  rows  or  columns, 
int  nRow()  const  {  return  nRow_;  } 
int  nCol()  const  {  return  nCol_;  } 

//  Parenthesis  operator  function. 

//  Allows  access  to  values  of  Matrix  via  (i,j)  pair. 

//  Example:  a(l,l)  =  2*b(2,3); 

//  If  column  is  unspecified,  take  as  1 . 
double&  operator))  (int  i,  int  j  =  1)  { 
assert(i  >  0  &&  i  <=  nRow_);  //  Bounds  checking  for  rows 
assert(j  >  0  &&j  <=  nCol  );  //  Bounds  checking  for  columns 

return  data_[  nCol_*(i-l)  +  (j-1)  ];  //  Access  appropriate  value 

} 


//  Parenthesis  operator  function  (const  version), 
const  double&  operator))  (int  i,  int  j  =  1)  const  ) 
assert(i  >  0  &&  i  <=  nRow_);  //  Bounds  checking  for  rows 
assert(j  >  0  &&j  <=  nCol  );  //  Bounds  checking  for  columns 

return  data_[  nCol_*(i-l)  +  (j-1)  ];  //  Access  appropriate  value 

} 


//  Set  function.  Sets  all  elements  of  a  matrix  to  a  given  value, 
void  set(double  value)  { 
int  i,  iData  =  nRow_*nCol_; 
for)  i=0;  i<iData;  i++  ) 
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data_[i]  =  value; 

} 


static  double  inv(Matrix  A,  Matrix&  Ainv)  { 

//  Compute  inverse  of  matrix 
//  Input 

//  A  -  Matrix  A  (N  by  N) 

//  Outputs 

//  Ainv  -  Inverse  of  matrix  A  (N  by  N) 
int  N  =  A.nRow(); 
assert(  N  ==  A.nCol() ); 

Ainv  =  A;  //  Copy  matrix  to  ensure  Ainv  is  same  size 
int  i,  j,  k; 

Matrix  scale(N),  b(N,N);  //  Scale  factor  and  work  array 
int  *index;  index  =  new  int  [N+l]; 

//*  Matrix  b  is  initialized  to  the  identity  matrix 
b.set(O.O); 

for(  i=l;  i<=N;  i++  ) 

b(i,i)=  10; 

//*  Set  scale  factor,  scale(i)  =  max(  |a(i,j)|  ),  for  each  row 
for(  i=l;  i<=N;  i++  )  { 

index[i]  =  i;  //  Initialize  row  index  list 

double  scalemax  =  0.; 

for(j=l;j<=N;j++) 

scalemax  =  (scalemax  >  fabs(A(i,j)))  ?  scalemax  :  fabs(A(i,j)); 
scale(i)  =  scalemax; 

} 


//*  Loop  over  rows  k  =  1, (N-l) 

int  signDet  =  1; 

for(  k=l;  k<=N-l;  k++  )  { 

//*  Select  pivot  row  from  max(  |a(j,k)/s(j)|  ) 
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double  ratiomax  =  0.0; 

int  jPivot  =  k; 

for(  i=k;  i<=N;  i++  )  { 

double  ratio  =  fabs(A(index[i],k))/scale(index[i]); 
if(  ratio  >  ratiomax  )  { 
jPivot=i; 
ratiomax  =  ratio; 

} 


//*  Perform  pivoting  using  row  index  list 
int  indexJ  =  index[k]; 
if(  jPivot  !=  k  )  {  //  Pivot 

indexJ  =  index[jPivot]; 

index[jPivot]  =  index[k];  //  Swap  index  jPivot  and  k 
index[k]  =  indexJ; 

signDet  *=  -1;  //  Flip  sign  of  determinant 


//*  Perform  forward  elimination 
for(  i=k+l;  i<=N;  i++  )  { 
double  coeff  =  A(index[i],k)/A(indexJ,k); 
for(j=k+l;j<=N;j++) 

A(index[i],j)  -=  coeff*  A(indexJ,j); 
A(index[i],k)  =  coeff; 
for(j=l;j<=N;j++) 

b(index[i],j)  -=  A(index[i],k)*b(indexJ,j); 


//*  Compute  determinant  as  product  of  diagonal  elements 
double  determ  =  signDet;  //  Sign  of  determinant 
for(  i=l;  i<=N;  i++  ) 

determ  *=  A(index[i],i); 

//*  Perform  backsubstitution 
for(  k=l;  k<=N;  k++  )  { 

Ainv(N,k)  =  b(index[N],k)/A(index[N],N); 
for(  i=N-l;  i>=l;  i~)  { 
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double  sum  =  b(index[i],k); 
for(j=i+l;j<=N;j++) 

sum  -=  A(index[i],j)*Ainv(j,k); 
Ainv(i,k)  =  sum/A(index[i],i); 

} 


delete  []  index;  //  Release  allocated  memory 
return(  determ ); 

} 


private: 

//  Matrix  data. 

int  nRow_,  nCol_;  //  Number  of  rows,  columns 

double*  data_;  //  Pointer  used  to  allocate  memory  for  data. 

//  Private  copy  function. 

//  Copies  values  from  one  Matrix  object  to  another. 

void  copy(const  Matrix&  mat)  { 
nRow_  =  mat.nRow_; 
nCol_  =  mat.nCol_; 
int  i,  iData  =  nRow_*nCol_; 
data_  =  new  double  [iData]; 
for(i  =  0;  i<iData;  i++  ) 
data_[i]  =  mat.data_[i]; 

} 


};  //  Class  Matrix 
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APPENDIX  C 


MATLAB  CODES  FOR  DOP  CALCULATION 


%  Code  to  obtain  DOP  (ENU)  from  Almanac  Data 
%  By  Yuen  Ming  Fatt 
%  Last  updated  on  27  Feb  2009 

9~'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 
o 

-k  ~k  -k  ~k  -k  -k  -k  ~k  -k  -k  -k  ~k  -k  -k  -k  -k  -k  -k  -k  -k  -k  -k  -k  ~k  -k  -k  -k  -k  -k  -k  -k  ~k 


•k'k-k'k-k'k-k'k-k'k-k'k-k'k-k'k-k'k-k'k-k'k-k'k-k'k-k'k'k'k'k'k'k 
: k'k'k'k'k'k-k'k'k'k'k'k’k'k-k'k'k'k'k'k'k'k-k'k-k'k'k'k'k'k-k'k'k 


clear  all 
clc 

format  long  g 


%********************************************************************** 
%Convert  Target's  Latitude,  Longitude  and  Altitude  to  ECEF  Coordinates 


%Target ' s  Latitude,  Longitude  and  Altitude  Input 
lat_deg  =  0;  %Latitude  (degree)  (user  input)##### 
lon_deg  =  90;  %Longitude  (degree)  (user  input)##### 
alt  =  0;  %Altitude  (meter)  (user  input)##### 
sealevel  =  0; 

lat  =  lat_deg . *pi . /I 80 ;  %Latitude  (rad) 

Ion  =  Ion  deg . *pi . /I 80 ;  %Longitude  (rad) 

Q, _ 

O 

%WGS84  ellipsoid  constants 
a  =  6378137; 

es  =  8.1819190842622e-2; 

o, _ 

o 


%Intermediate  calculation 

N  =  a. /sqrt (1-es . A2 . *sin (lat) . A2) 

%Results 

xTgt  =  (N+alt) . *cos (lat) . *cos (Ion 
yTgt  =  (N+alt) . *cos (lat) . *sin (Ion 
zTgt  =  (  ( 1-es . A2 )  . *N  +  alt).*sin( 
(meter) 


;  %Prime  vertical  radius  of  curvature 

) ;  %Target  x  ECEF  coordinate  (meter) 

) ;  %Target  y  ECEF  coordinate  (meter) 
lat) ;  %Target  z  ECEF  coordinate 


xLocalRef  =  (N+sealevel 
coordinate  (meter) 
yLocalRef  =  (N+sealevel 
coordinate  (meter) 
zLocalRef  =  ((l-es.A2). 
ECEF  coordinate  (meter) 
%End  of  conversion 


) . *cos (lat) . * 

cos (Ion) ; 

%ENU  Local 

ref  pt 

X 

ECEF 

) . *cos (lat) . * 

sin (Ion) ; 

%ENU  Local 

ref  pt 

y 

ECEF 

*N  +  sealevel 

) . *sin (lat 

) ;  %ENU  Local  ref 

pt 

z 

^^it^^ititit^itit^it'k'k'k'k'k'k'k'k'k-k'k'k'k'k'k'k'k'k'k'k-k'k'k'k-k'k'k'k-k'k'k'k-k'k'k'k-k'k'k'k-k'k'k'k-k'k'k'k'k'k'k'k-k'k'k'k'k'k 

o 

%Convert  Almanac  Data  to  Satellite  Position  in  ECEF  Coordinates  and 
check  the  Line  of  Sight  to  Target 

★  -k  -k  -k  -k  -k  -k  -k  ~k  -k  -k  -k  ~k  -k  -k  -k  ★  -k  -k  -k  ~k  -k  -k  -k  -k  -k  -k  -k  ~k  -k  -k  -k  ~k  -k  -k  -k  -k  -k  -k  -k  -k  -k  -k  -k  ★  -k  -k  -k  ~k  -k  -k  -k  ★  -k  -k  -k  ~k  -k  -k  -k  -k  -k  -k  -k  ~k  -k  -k  -k  -k  -k 

%Constants 

io  =  0.3.*pi;  %Inclination  angle  @  ref.  time  (rad) 

mju  =  3.986005el4;  %WGS  84  value  of  the  Earth's  universal  gravitational 
parameter  for  GPS  user  (metersA3/secA2 ) 
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OMEGAdote  =  7 . 2 92 1 1514 67e-5 ;  %WGS  84  value  of  the  Earth's  rotation  rate 
(rad/sec) 

o, _ 

o 

%Almanac  Data  from  Satellite 

fid  =  f open (' current . al3 ') ;  %Open  source  file  "current . al3" 

%Common  Data 

NumSV  =  fscanf(fid,  '%d');  %Number  of  Satellites 
name  =  fgetl (fid) ; 

[data,  count]  =  f scanf ( f id,  ' %f  ' )  ; 
fclose (fid) ; 

wn  =  data(l);  %GPS  week  no. 

toa  =  data (2);  %Time  of  Applicability  of  Almanac (sec)  (range:  0  to 
604,784) 


o, _ 

o 

%Input  time  from  user##### 

y  =  2009;  %Year 

m  =  2;  %Month 

d  =  3;  %Date 

h  =  0;  %Hours 

mi  =  0;  %Minutes 

sec  =  0;  %Seconds 

timezone  =  0;  %Timezone  (Eastern  Standard  Time  (North  America)  =  -5hr) 
summertime  =  0;  %To  account  for  daylight  saving.  If  summer,  1  =  Yes,  0 
=  No 


o. _ 

o 

[gps  week,  sec  of  week]  =  ymdhms2gps (y,  m,  d,  h,  mi,  sec,  timezone, 
summertime) ; 

Total  weeks  =  gps  week;  %Total  number  of  weeks  since  6  Jan  1980 
while  gps  week  >=  1024 

gps  week  =  gps  week  -  1024; 

end 

tk  =  (gps_week  -  wn) *604800  +  (sec_of_week  -  toa);  %Time  since  toa (sec) 
(range:  -302,400  to  302,400) 
if  gps  week  <  wn 

fprintf (' Almanac  file  used  is  incorrect.  Please  use  almanac  file 
for  week  %4.0f',gps  week) 
elseif  gps  week  >  wn 

fprintf (' Almanac  file  is  outdated.  Please  use  almanac  file  for  week 
%4 . Of ' , gps_week) 
end 


o, _ 

o 

%Satellite  Specific  Data 
num  =  0 ; 

SVcount  =  1; 

while  SVcount  <=  NumSV, 

PRN  =  data(num+3);  %PRN  number 

SVN  =  data(num+4);  %Satllite  number 

URA  =  data(num+5);  %  Average  URA  number 

ec  =  data(num+6);  %Eccentricity  (dimensionless)  (range:  0-0.03) 
del  ik  =  data(num+7) .*pi;  %Inclination  correction  (rad) 

OMEGAdot  =  data(num+8) .*pi;  %Rate  of  right  ascension  (rad/sec) 
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sqrtA  =  data(num+9);  %Sqr  root  semi-major  axis  (mAl/2) 

OMEGAo  =  data(num+10) . *pi;  %Right  ascension  @  ref.  time  (rad) 
omega  =  data(num+ll) . *pi;  %Argument  of  perigee  (rad) 

Mo  =  data(num+12) . *pi;  %Mean  Anomaly  @  ref.  time  (rad) 

AfO  =  data (num+13 ) ;  %Clock  offset  (sec) 

Afl  =  data (num+14 ) ;  %Clock  drift  (sec/sec) 

Health  =  data (num+15) ;  %Satellite  Health;  0=healthy 
num  =  num+14; 

%End  of  data  extraction 

O, _ 

o 

%Calculations 

A  =  sqrtA. A2;  %Orbit  semi-major  axis  (meter) 
n  =  sqrt (mj u . / (A . A3 ) ) ;  %Computed  mean  motion  (rad/sec) 

Mk  =  Mo+tk.*n;  %Mean  anomaly  (rad) 

%Start  values  for  iterative  solution  of  Kepler  eq. 

Ek  =  Mk; 

Eold=0 ; 

while  abs (Ek-Eold) >=1 . 0e-10 
Eold  =  Ek; 

Ek  =  Mk+ec . *sin (Ek) ; 

end 

%End  of  iteration 

vk  =  atan2 ( ( sqrt (l-ec.A2)  . *sin (Ek) )  ./ (l-ec.*cos(Ek)),  (cos (Ek) - 
ec)  . / (1-ec . *cos (Ek) ) ) ;  %True  anomaly  (rad) 

Ek  =  acos  (  (ec+cos (vk) )  . / ( 1+ec . *cos (vk) ) ) ; 

uk  =  omega+vk;  %Argument  of  latitude  (rad) 

rk  =  A. * ( 1-ec . *cos (Ek) ) ;  %Corrected  radius  (meter) 

ik  =  io+del  ik;  %Corrected  inclination  (rad) 

xkl  =  rk.*cos(uk);  %x  position  in  orbital  plane  (meter) 

ykl  =  rk.*sin(uk);  %y  position  in  orbital  plane  (meter) 

OMEGAk  =  OMEGAol (OMEGAdot-OMEGAdote) . *tk-OMEGAdote . * (toa) ; 
%Corrected  longitude  of  ascending  node  (rad) 

o, _ 

o 

%Calculations  for  ECEF  coordinates 

xk(SVcount)  =  xkl . *cos (OMEGAk) -ykl . *cos (ik) . *sin (OMEGAk) ; 
%Satellite  x  ECEF  coordinate  (meter) 

yk(SVcount)  =  xkl . *sin (OMEGAk) +ykl . *cos (ik) . *cos (OMEGAk) ; 
%Satellite  y  ECEF  coordinate  (meter) 

zk(SVcount)  =  ykl . *sin (ik) ;  %Satellite  z  ECEF  coordinate  (meter) 
%End  of  ECEF  coordinates  conversion  from  Almanac  Data 

O, _ 

o 

%Covert  ECEF  coordinates  to  East-North-Up  Coordinates 
East  =  -sin  ( Ion)  . * (xk-xLocalRef )  +  cos  (Ion)  . * (yk-yLocalRef ) ; 

North  =  -sin (lat) . *cos (Ion) . * (xk-xLocalRef )  - 
sin (lat) . *sin (Ion) . * (yk-yLocalRef )  +  cos (lat) . * (zk-zLocalRef ) ; 

Up  =  cos (lat)  . *cos (Ion)  . * (xk-xLocalRef )  +  cos (lat)  . *sin (Ion)  . * (yk- 
yLocalRef)  +  sin (lat) . * (zk-zLocalRef ) ; 

o, _ 

o 

%Algorithm  to  determine  Line  of  Sight  between  Target  and  Satellite 
Obstruction  =  10;  %Obstruction  to  the  Field  of  View  (deg)  (user 
input) ##### 

mag  Tgt  =  sqrt (xTgt . A2+yTgt . A2  +  zTgt . A2 )  ;  %Distance  of  Target  from 
Earth  center  (meter) 

xTgttoSV (SVcount)  =  xk(SVcount)  -  xTgt; 
yTgttoSV (SVcount)  =  yk (SVcount)  -  yTgt; 
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zTgttoSV (SVcount)  =  zk(SVcount)  -  zTgt; 
mag_TgttoSV (SVcount)  = 

sqrt ( (xTgttoSV (SVcount) ) . A2  + (yTgttoSV (SVcount) ) . A2+ (zTgttoSV (SVcount) ) . 
A2);  %Distance  from  Target  to  Satellite 

AngleFromTgt (SVcount)  =  acos (( (xTgttoSV (SVcount) . *xTgt)  + 

(yTgttoSV (SVcount) . *yTgt)  + 

(zTgttoSV (SVcount) . *zTgt) ) . / (mag_TgttoSV (SVcount) . *mag_Tgt) ) ; 


%if  mag_SVproj (SVcount) >mag  Tgt  &&  AngleTOS (SVcount) < (pi/2)  && 

AngleFromTgt (SVcount) < (pi/2- (Obstruction*pi/180) )  &&  Health  ==0 

if  AngleFromTgt (SVcount) < (pi/2- (Obstruction*pi/180) )  &&  Health 

Los  =  1;  %There  is  Line  of  Sight 

else 

Los  =  0;  %There  is  NO  Line  of  Sight 

end 

LOS (SVcount)  =  Los; 

SVcount  =  SVcount+1; 

end 


==0 


o 

%Assign  coordinates  to  Target  and  valid  Satellites  (i.e.  those  with  LOS 
with  Target) 

o 

num3  =  1 ; 
num4  =  1 ; 

while  num3  <=  NumSV 
if  L0S(num3)  ==  1 

SV(num4,l:3)  =  [East(num3),  North (num3),  Up(num3)];  %Assigning 
coordinates  to  respective  valid  satellites 
num4  =  num4+l; 
end 

num3  =  num3+l; 

end 

NumValidSV  =  num4-l; 


9~~k  -k  ~k  -k  ~k  -k  ~k  -k  ~k  -k  ~k  -k  ~k  -k  ~k  -k  ~k  -k  ~k  -k  ~k  -k  ~k  -k  ~k  -k  ~k  -k  ~k  -k  ~k  -k  ~k  -k  ~k  -k  ~k  -k  ~k  -k  ~k  -k  ~k  -k  ~k  -k  ~k  -k  ~k  -k  ~k  -k  ~k  -k  ~k  -k  ~k  -k  ~k  -k  ~k  -k  ~k  -k  ~k  -k  ~k  -k  ~k  -k 
o 

%Calculate  PDOP  with  the  valid  satellites  i.e.  those  within  the  LOS  of 
Tgt 

~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k 

%Pseudo-Range  and  Directional  Derivative  Loop 
for  num4  =  1: NumValidSV 

%Calculate  pseudo-ranges  from  reciever  position  to  other  vehicles 
r(num4)  =  sqrt ( ( SV (num4 , 1 ) ) A2  +  ( SV (num4 , 2 ) ) A2  +  (SV(num4,3)- 
alt) A2 )  ; 

%Calculate  directional  derivatives  for  X,Y,Z,  and  Time 
Dx(num4)  =  ( SV (num4 , 1 ) -0 ) /r (num4 ) ;  %x-coordinates  of  Tgt  in  ENU 
frame  is  zero 

Dy(num4)  =  ( SV (num4 , 2 ) -0 ) /r (num4 ) ;  %y-coordinates  of  Tgt  in  ENU 
frame  is  zero 

Dz(num4)  =  ( SV (num4 , 3 ) -alt ) /r (num4 ) ;  %z-coordinates  of  Tgt  in  ENU 
frame  is  the  altitude 
Dt (num4)  =  -1; 
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end 


%Produce  the  Covariance  Matrix  from  the  Directional  Derivatives 
Alp  =  zeros (NumValidSV, 4 ) ; 
for  num5  =  l:NumValidSV 

Alp(num5,l)  =  Dx(num5); 

Alp(num5,2)  =  Dy(num5); 

Alp(num5,3)  =  Dz(num5); 

Alp(num5,4)  =  Dt(num5); 

end 

Brv  =  transpose (Alp) ; 

Chi  =  Brv*Alp; 

Dlt  =  inv (Chi ) ; 


%  Calculate  DOPs  from  the  diagonal  elements  of  the  Covariance  Matrix 
GDOP  =  sqrt (Dlt (1, 1)  +  Dlt(2,2)  +  Dlt(3,3)  +  Dlt(4,4)); 

PDOP  =  sqrt (Dlt (1, 1)  +  Dlt(2,2)  +  Dlt(3,3)); 

HDOP  =  sqrt (Dlt (1, 1)  +  Dlt(2,2)); 

TDOP  =  sqrt (Dlt  (4, 4) )  ; 

VDOP  =  sqrt (Dlt (3, 3) ) ; 

YDOP  =  sqrt (Dlt  (2,2)); 

XDOP  =  sqrt (Dlt  (1, 1) )  ; 


Results  ( 1 )  = 
Results  (2 )  = 
Results  ( 3 )  = 
Results  ( 4 )  = 
Results  ( 5 )  = 
Results  ( 6 )  = 
Results  ( 7 )  = 
Results  ( 8 )  = 

%End  of  Code 


NumValidSV; 

GDOP; 

PDOP; 

HDOP; 

TDOP; 

VDOP; 

YDOP; 

XDOP; 

:) 


158 


THIS  PAGE  INTENTIONALLY  LEFT  BLANK 


159 


VDOP  .  VDOP  VDOP 


APPENDIX  D.  DOP  COMPARISON  BETWEEN  MATLAB  & 

TRIMBLE 
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dOQA  dOQA  dOQA 
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B 


TDOP 


TDOP  at  North  &  South  Pole 
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150  2-87  253  1.15  135  225  9  2370  2534  1.149  1348  2.258  0319  0305  9  0.00  0.14  -0.11  -0.15  037 


Lat<-f/-90|  Long{+/-180|  z{  \z  z{  Difference  m 

+North/-South|  {+East/-West|  GPOP  PDOP  HDOP  TDOP  VDOP  Visible  SV  GDOP  PDOP  HDOP  TDOP  VDOP  YD  OP  XDOP  VgibteSV  GDOP  PDOP  HDOP  TDOP  VDOP  Visible  SV 
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165  156  141  056  0.67  111  9  1560  1409  0562  0.670  1.115  0.622  0597  9  0.01  -0.06  0 24  -0.04  041 


Lat{4-/-90)  Long{+/-180)  No. of  Difference  in 

4North/-South|  {4  East/ -West)  6D0P  PDOP  HDOP  TDOP  VDOP  Visible  SV  GDOP  PDOP  HDOP  TDOP  VDOP  YD  OP  XDOP  Visible  SV  GDOP  PDOP  HDOP  TDOP  VDOP  Visible  SV 
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180  1.98  1.75  0.94  C.93  147  10  1.955  1.726  0.938  0.918  1449  0.773  0531  10  1-1-29  -138  -024-133 


Lat{+/-90)  long{+/-180|  \:  No. of  Difference  in 

4North/-South|  {+East/-West|  GDOP  PDOP  HDOP  TDOP  VDOP  Visible  SV  GDOP  PDOP  HDOP  TDOP  VDOP  YD  OP  XDOP  Visible  SV  GDOP  PDOP  HDOP  TDOP  VDOP  Visible  SV 
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lac  153  1.63  0.7a  054  143  11  1513  1.614  0.779  0526  1413  0572  0529  11  -0.95  -1.00  -O.Ca  -1.70  -1.19 


Lat{+/-90|  Long{+/-180|  No.  a#  \z  z{  Difference  m 

North/»Sauth|  (+EasV-West|  GDOP  PDOP  HDOP  TDOP  VDQP  Visible  SV  GPOP  PDOP  HDOP  TDOP  VDOP  YD  OP  XDOP  Visible  SV  GDOP  PDOP  HDOP  TDOP  VDOP  Visible  SV 
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155  2.01  1.73  0.73  0.93  1.6  11  2.034  1,802  0.785  0.944  1.622  0-521  0-583 


Lat{+/-9C)  Long(+/-180]  \:  :*  No. of  Difference  in 

+North/-South|  {+EasV*West|  GDOP  PDOP  HDOP  TDOP  VDOP  Visible  SV  GDOP  PDOP  HDOP  TDOP  VDOP  YDOP  XDOP  Visible  SV  GDOP  PDOP  HDOP  TDOP  VDOP  Visible  SV 
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150  1.77  153  0.78  0.78  138  11  1.785  1.601  0.780  0.791  1398  0580  0521 


No. of  No. of  I  Difference  in 

PDOP  HDOP  TDOP  VDOP  Visible  SV  GDOP  PDOP  HDOP  TDOP  VDOP  YDOP  XDOP  Visible  SV  IgDOP  PDOP  HDOP  TDOP  VDOP  Visible  SV 
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165  1.94  1.72  1.05  0.9  136  8  1.932  1.712  1.052  0.894  1351  0.713  0.773  8  -043  -045  0.15  -0.66  -0.65 


Lat{+/-90|  Long{+/-180|  \:  =?  Difference  m 

North/-5outh|  {+East/-West|  GDOP  PDOP  HDOP  TDOP  VDOP  Visible  SV  GDOP  PDOP  HDOP  TDOP  VDOP  YDOP  XDOP  Visible  SV  k;DOP  PDOP  HDOP  TDOP  VDOP  Visible  SV 
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ISO  233  2.03  l.oa  1.15  1.71  3  2387  1.993  1.082  1.122  1.674  0353  0.665  3  1-133  -131  0.16  -241  -2.10 


Lat{+/-90|  Longj-fZ-iaOj  Na  of  \:  :*  Difference  m 

North/-5outh|  (»East/-West|  GDOP  PDOP  HDOP  TDOP  VDOP  Visible  SV  GDOP  PDOP  HDOP  TDOP  VDOP  YD  OP  XDOP  Visible  SV  GDOP  PDOP  HDOP  TDOP  VDOP  Visible  SV 
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ISO  2.13  153  053  1  1.66  9  2.093  1551  0577  0.933  1.630  0.642  0597  9  1-152  -156  -059  -155  -152 


Lat{+/-90)  Long{+/-180|  No. of  No. of  Difference  in 

+North/-South|  {+ East/- West)  GDOP  PDOP  HDOP  7DOP  VDOP  VisWeSV  GDOP  PDOP  HDOP  TDOP  VDOP  YDOP  XDOP  Visible  SV  GDOP  PDOP  HDOP  TDOP  VDOP  Visible  SV 
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18C  1.93  174  0.33  053  153  9  1.901  1.716  0535  0517  1500  0561  0.613  9  1-152  -156  059  -1.63  -1.99 


Lat{4/-90)  Long{4/-13C|  No. of  No. of  Difference  in 

+North/-5outh|  J 4 East/ -West)  IgDOP  PDOP  HDOP  TDOP  VDOP  Visible  S^’  GDOP  PDOP  HDOP  TDOP  VDOP  YDOP  XDOP  Visible  SV  GDOP  PDOP  HDOP  TDOP  VDOP  Visible  SV 
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PDOP  ,  PDOP  PDOP 


D 


PDOP 


PDOP  at  North  &  South  Pole 


i 

o 

2.5 

,  □ 

□ 

2 

1.5 

1 

0.5 

— 0-4  — , —  — , 

-90  0  90 


PDOP  vs  Longitude  (deg)  for 
Latitude  -60  deg 


3 


O.  1 


0.5  A 
0 

-180-150-120  -90  -60  -30  0  30  60  90  120  150  180 


Latitude  (deg)  (+North/  -South)  Longitude  (deg)  (+East/  -West) 


O  Trimble  □  Visual  C++ 


o  Trimble  □  Visual  C++ 


PDOP  vs  Longitude  (deg)  for 
Latitude  -30  deg 


3  j 

G 

2.5 

3  O 

a 

0  □  □  2 
□  1.5 

□  n 

□  Q  H  E 

1 

0.5 

1 - I - 1 - 1 - I - T0— 

1-150-120  -90  -60  -30  0  30  60  90  120  150  180 


Longitude  (deg)  (+East/  -West) 

o  Trimble  □  Visual  C++ 


PDOP  vs  Longitude  (deg)  for 
Latitude  0  deg 

2.5  n 


0.5 

-T&- 


-180-150-120  -90  -60  -30  0  30  60  90  120  150  180 

Longitude  (deg)  (+East/  -West) 

o  Trimble  □  Visual  C++ 


PDOP  vs  Longitude  (deg)  for  PDOP  vs  Longitude  (deg)  for 

Latitude  30  deg  Latitude  60  deg 


i - 1 - 1 - 1 - 1 - r-e — I 1 - 1 - 1 - 1 - 1 - 1  i - 1 - 1 - 1 - 1 - r0  I - 1 - 1 - 1 - 1 - 1 — , 

-180-150-120  -90  -60  -30  0  30  60  90  120  150  180  -180-150-120  -90  -60  -30  0  30  60  90  120  150  180 

Longitude  (deg)  (+East/  -West)  Longitude  (deg)  (+East/  -West) 

o  Trimble  □  Visual  C++  o  Trimble  □  Visual  C++ 


185 


GDOP  .  GDOP  GDOP 
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•150  133  1.5a  032  C.33  147  9  1374  1.584  0318  0322  1471  059  0558  9  -032  024  -024  -0.95  0.07 
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180  11.83  1.63  0.78  034  143  11  |l313  1.614  0.78  0326  1413  0573  0529  11  |-0.93  -0.98  0.00  -1.57  -1.19 
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APPENDIX  F.  COMPARISON  OF  VISIBLE  SATELLITES 
BETWEEN  MATLAB  &  TRIMBLE 
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Angle  between  GPS  Receiver  &  Satellite  (degrees)  wrt  vertical 
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APPENDIX  G.  COMPARISON  OF  DOP  VALUES  FROM 
OUTDATED  ALMANAC  DATA 
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APPENDIX  H.  AVERAGE  DOP  VALUES 


Date 

Jul  28,  2008 

Jul  31,  2008 

Aug  11,  2008  Aug  27,  2008 

Dec  3,  2008 

Feb  3,  2009 

Global  Average 

XDOP 

0.516 

0.516 

0.523 

0.529 

0.514 

0.515 

0.519 

YDOP 

0.55 

0.55 

0.558 

0.565 

0.548 

0.549 

0.553 

Average  DOP 

VDOP 

1.202 

1.201 

1.221 

1.242 

1.194 

1.197 

1.210 

for 

Obstruction 

TDOP 

0.649 

0.649 

0.66 

0.672 

0.644 

0.646 

0.653 

Angle  =  0° 

HDOP 

0.758 

0.758 

0.768 

0.777 

0.755 

0.757 

0.762 

PDOP 

1.426 

1.426 

1.447 

1.47 

1.417 

1.421 

1.435 

GDOP 

1.567 

1.567 

1.591 

1.617 

1.557 

1.562 

1.577 

XDOP 

0.637 

0.636 

0.636 

0.635 

0.632 

0.633 

0.635 

YDOP 

0.715 

0.714 

0.713 

0.711 

0.71 

0.711 

0.712 

Average  DOP 

VDOP 

1.772 

1.77 

1.767 

1.762 

1.758 

1.759 

1.765 

for 

Obstruction 

TDOP 

1.08 

1.079 

1.076 

1.073 

1.071 

1.071 

1.075 

Angle  =  10° 

HDOP 

0.966 

0.966 

0.964 

0.962 

0.959 

0.96 

0.963 

PDOP 

2.029 

2.027 

2.023 

2.017 

2.013 

2.014 

2.021 

GDOP 

2.3 

2.298 

2.293 

2.287 

2.282 

2.283 

2.291 

XDOP 

0.732 

0.735 

0.743 

0.756 

0.725 

0.721 

0.735 

YDOP 

0.885 

0.89 

0.894 

0.918 

0.865 

0.853 

0.884 

Average  DOP 

VDOP 

2.361 

2.393 

2.417 

2.512 

2.319 

2.262 

2.377 

for 

Obstruction 

TDOP 

1.542 

1.564 

1.581 

1.65 

1.508 

1.47 

1.553 

Angle  =  15° 

HDOP 

1.169 

1.175 

1.182 

1.209 

1.145 

1.132 

1.169 

PDOP 

2.654 

2.686 

2.712 

2.811 

2.605 

2.546 

2.669 

GDOP 

3.072 

3.111 

3.143 

3.263 

3.014 

2.962 

3.094 
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APPENDIX  I.  COMPARISON  BETWEEN  VDOP  &  HDOP  FOR 

OBSTRUCTION  ANGLE  =  10° 


Lat  (+/-90) 

From  Matlab 

Long  (+/-180)  (+East/-West) 

(+North/-South) 

Program 

-165 

-150 

-135 

-120 

-105 

-90 

-75 

-60 

-45 

-30 

-15 

0 

HDOP 

1.0625 

0.9999 

0.9639 

0.7899 

0.7551 

0.7345 

0.7625 

0.7765 

0.8121 

0.8432 

0.9372 

0.9481 

-75 

VDOP 

2.2301 

2.1614 

2.1693 

1.6064 

1.6105 

1.6104 

1.8262 

1.815 

1.8024 

1.8412 

2.081 

2.0896 

VDOP/HDOP 

2.10 

2.16 

2.25 

2.03 

2.13 

2.19 

2.40 

2.34 

2.22 

2.18 

2.22 

2.20 

HDOP 

0.9585 

0.8182 

0.7652 

0.88 

0.8897 

0.9376 

0.7838 

0.8109 

0.8257 

0.7788 

1.0127 

1.1056 

-60 

VDOP 

1.5495 

1.4715 

1.399 

1.7528 

1.7527 

1.8916 

1.3482 

1.3835 

1.3694 

1.1682 

1.2553 

1.3884 

VDOP/HDOP 

1.62 

1.80 

1.83 

1.99 

1.97 

2.02 

1.72 

1.71 

1.66 

1.50 

1.24 

1.26 

HDOP 

0.7845 

1.0444 

0.8567 

0.8024 

0.9586 

1.1863 

1.0722 

0.9991 

0.9872 

0.9932 

1.2276 

1.1133 

-45 

VDOP 

1.226 

1.9465 

1.719 

1.5231 

1.6833 

1.5174 

1.6146 

1.339 

1.3713 

1.3789 

1.877 

1.6956 

VDOP/HDOP 

1.56 

1.86 

2.01 

1.90 

1.76 

1.28 

1.51 

1.34 

1.39 

1.39 

1.53 

1.52 

HDOP 

0.8891 

1.1133 

0.94 

0.8012 

0.9545 

1.126 

1.0938 

1.0556 

1.0681 

1.1437 

0.8942 

1.0076 

-30 

VDOP 

1.5539 

1.6167 

1.5316 

1.3761 

1.5711 

1.5216 

1.3059 

1.5371 

1.5207 

1.8955 

1.5889 

2.5775 

VDOP/HDOP 

1.75 

1.45 

1.63 

1.72 

1.65 

1.35 

1.19 

1.46 

1.42 

1.66 

1.78 

2.56 

HDOP 

0.8969 

0.9465 

0.8757 

1.0302 

0.9421 

0.7521 

1.0238 

0.9769 

1.0051 

0.9153 

0.7738 

0.9584 

-15 

VDOP 

1.6166 

1.5739 

1.3769 

1.6733 

1.5746 

1.3061 

1.7545 

1.5795 

1.5527 

1.4673 

1.321 

2.6471 

VDOP/HDOP 

1.80 

1.66 

1.57 

1.62 

1.67 

1.74 

1.71 

1.62 

1.54 

1.60 

1.71 

2.76 

HDOP 

0.8261 

0.9784 

0.8688 

0.9202 

0.7573 

0.758 

0.8037 

0.9542 

0.7587 

0.8196 

0.8464 

0.8627 

0 

VDOP 

1.2304 

1.5315 

1.4217 

1.845 

1.4634 

1.479 

1.6135 

1.6283 

1.0743 

1.3739 

1.3908 

1.3948 

VDOP/HDOP 

1.49 

1.57 

1.64 

2.00 

1.93 

1.95 

2.01 

1.71 

1.42 

1.68 

1.64 

1.62 

HDOP 

1.0409 

0.8134 

0.9886 

0.9577 

0.981 

0.7759 

0.8367 

0.7649 

0.9597 

0.7841 

0.9984 

0.8391 

15 

VDOP 

1.4624 

1.0984 

1.6328 

2.3865 

2.4136 

1.5255 

2.0938 

1.3559 

1.7049 

1.2506 

1.756 

1.2575 

VDOP/HDOP 

1.40 

1.35 

1.65 

2.49 

2.46 

1.97 

2.50 

1.77 

1.78 

1.59 

1.76 

1.50 

HDOP 

1.0529 

1.1707 

0.8603 

0.9373 

0.8881 

0.9799 

1.1539 

0.9204 

0.8318 

0.9465 

0.8983 

1.153 

30 

VDOP 

1.8536 

1.8373 

1.3208 

1.6216 

1.6382 

1.5464 

2.3453 

1.6779 

1.2521 

1.5605 

1.4153 

1.7551 

VDOP/HDOP 

1.76 

1.57 

1.54 

1.73 

1.84 

1.58 

2.03 

1.82 

1.51 

1.65 

1.58 

1.52 

HDOP 

0.8857 

1.1367 

0.9718 

1.0757 

1.1747 

1.211 

1.8173 

1.1202 

1.1452 

1.2925 

0.9871 

0.8529 

45 

VDOP 

1.6452 

2.7577 

1.5247 

1.5603 

1.5268 

1.3642 

2.6174 

1.5663 

1.5405 

2.3618 

1.9303 

1.5966 

VDOP/HDOP 

1.86 

2.43 

1.57 

1.45 

1.30 

1.13 

1.44 

1.40 

1.35 

1.83 

1.96 

1.87 

HDOP 

0.8329 

0.7672 

0.7985 

1.0052 

1.81 

1.9523 

1.503 

1.5491 

0.955 

0.7997 

0.919 

0.8921 

60 

VDOP 

1.6323 

1.5154 

1.4979 

1.7956 

3.2854 

3.192 

1.258 

1.6089 

1.3477 

1.3313 

1.7181 

1.9745 

VDOP/HDOP 

1.96 

1.98 

1.88 

1.79 

1.82 

1.63 

0.84  | 

1.04 

1.41 

1.66 

1.87 

2.21 

HDOP 

0.7976 

0.8176 

1.0087 

1.0545 

1.1446 

1.1119 

0.9453 

0.9397 

0.9897 

0.9589 

0.8009 

0.7692 

75 

VDOP 

2.0214 

2.0239 

3.6001 

3.5553 

3.5077 

2.6718 

2.0048 

2.0016 

2.3267 

2.3437 

1.7415 

1.7577 

VDOP/HDOP 

2.53 

2.48 

3.57 

3.37 

3.06 

2.40 

2.12 

2.13 

2.35 

2.44 

2.17 

2.29 

HDOP 

0.7752 

90 

VDOP 

2.0945 

VDOP/HDOP 

2.70 

HDOP 

0.7591 

-90 

VDOP 

1.7431 

VDOP/HDOP 

2.30 
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Lat  (+/-90) 
(+North/-South) 

From  Matlab 
Program 

!  Long  (+/-180)  (+East/-West) 

-165 

-150 

-135 

-120 

-105 

-90 

-75 

-60 

-45 

-30 

-15 

0 

HDOP 

0.7944 

0.8576 

0.8083 

0.7838 

0.788 

0.8201 

0.875 

0.9452 

1.0232 

1.1487 

0.9071 

0.9576 

-75 

VDOP 

1.716 

2.1094 

2.1104 

2.106 

2.0962 

2.0816 

2.0632 

2.0425 

2.0214 

2.2583 

1.6538 

1.6753 

VDOP/HDOP 

2.16 

2.46 

2.61 

2.69 

2.66 

2.54 

2.36 

2.16 

1.98 

1.97 

1.82 

1.75 

HDOP 

1.1659 

1.1252 

1.0991 

1.0589 

0.805 

0.8183 

0.8726 

0.8966 

0.9209 

0.9718 

0.8621 

0.8341 

-60 

VDOP 

1.6762 

1.5156 

2.0696 

2.0705 

1.4313 

1.4146 

1.3852 

1.3301 

1.3101 

1.2843 

1.1146 

1.1732 

VDOP/HDOP 

1.44 

1.35 

1.88 

1.96 

1.78 

1.73 

1.59 

1.48 

1.42 

1.32 

1.29 

1.41 

HDOP 

0.8257 

0.8862 

0.8058 

0.8757 

1.0216 

1.0521 

1.0043 

0.982 

1.0391 

1.1071 

1.4274 

0.9377 

-45 

VDOP 

1.287 

1.357 

1.2419 

1.3952 

1.7959 

1.7727 

1.4092 

1.412 

1.3863 

1.369 

2.3315 

1.4488 

VDOP/HDOP 

1.56 

1.53 

1.54 

1.59 

1.76 

1.68 

1.40 

1.44 

1.33 

1.24 

1.63 

1.55 

HDOP 

0.8592 

0.9217 

0.9002 

0.8921 

0.9812 

1.0376 

0.9853 

0.8617 

0.9491 

1.1597 

0.9738 

0.7794 

-30 

VDOP 

1.7422 

1.3812 

1.4426 

1.767 

1.7247 

1.4541 

1.4776 

1.2248 

1.5522 

2.4277 

2.0797 

1.413 

VDOP/HDOP 

2.03 

1.50 

1.60 

1.98 

1.76 

1.40 

1.50 

1.42 

1.64 

2.09 

2.14 

1.81 

HDOP 

0.7592 

0.9618 

0.9231 

0.773 

0.9219 

0.933 

0.847 

0.9502 

0.8104 

0.7251 

0.7854 

1.0624 

-15 

VDOP 

1.3859 

1.896 

2.1601 

1.4265 

1.5445 

1.3396 

1.1934 

1.6841 

1.3686 

1.3572 

1.6216 

2.3856 

VDOP/HDOP 

1.83 

1.97 

2.34 

1.85 

1.68 

1.44 

1.41 

1.77 

1.69 

1.87 

2.06 

2.25 

HDOP 

0.7275 

0.8056 

0.8379 

0.7659 

0.9582 

0.8064 

0.8565 

0.7197 

0.766 

0.7801 

0.8587 

0.9118 

0 

VDOP 

1.2525 

1.5234 

1.6099 

1.3801 

1.6441 

1.1709 

1.6742 

1.4733 

1.4483 

1.3977 

1.5548 

1.4238 

VDOP/HDOP 

1.72 

1.89 

1.92 

1.80 

1.72 

1.45 

1.95 

2.05 

1.89 

1.79 

1.81 

1.56 

HDOP 

0.8667 

0.7395 

1.0045 

0.9866 

0.7814 

0.8438 

0.9131 

0.7395 

0.9047 

0.7803 

1.0516 

1.0147 

15 

VDOP 

1.4213 

1.1908 

1.7332 

1.7296 

1.0833 

1.3334 

2.445 

1.7692 

1.8635 

1.273 

1.3512 

1.5009 

VDOP/HDOP 

1.64 

1.61 

1.73 

1.75 

1.39 

1.58 

2.68 

2.39 

2.06 

1.63 

1.28 

1.48 

HDOP 

1.1499 

0.9721 

1.0491 

1.0291 

0.9251 

0.9652 

1.2426 

0.9292 

1.0271 

1.2077 

1.1479 

1.0817 

30 

VDOP 

1.5707 

1.266 

1.514 

1.4928 

1.284 

1.3687 

2.6182 

2.1491 

2.082 

1.6792 

1.6341 

1.6741 

VDOP/HDOP 

1.37 

1.30 

1.44 

1.45 

1.39 

1.42 

2.11 

2.31 

2.03 

1.39 

1.42 

1.55 

HDOP 

0.8955 

1.0091 

1.0674 

1.065 

0.9418 

1.4463 

1.4629 

1.311 

1.3456 

0.9346 

0.9773 

0.8766 

45 

VDOP 

1.6025 

1.3686 

1.8086 

1.7611 

1.2625 

1.7182 

2.0918 

2.1355 

2.0968 

1.1849 

1.5796 

1.6298 

VDOP/HDOP 

1.79 

1.36 

1.69 

1.65 

1.34 

1.19 

1.43 

1.63 

1.56 

1.27 

1.62 

1.86 

HDOP 

0.785 

0.9297 

0.8529 

0.9995 

1.0755 

1.2293 

0.9602 

0.8795 

0.8616 

0.838 

0.8327 

0.8349 

60 

VDOP 

1.7413 

2.2768 

1.6408 

1.8467 

1.793 

1.81 

1.3778 

1.3048 

1.3145 

1.3845 

1.4293 

1.4995 

VDOP/HDOP 

2.22 

2.45 

1.92 

1.85 

1.67 

1.47 

1.43 

1.48 

1.53 

1.65 

1.72 

1.80 

HDOP 

0.8595 

0.8575 

0.9034 

0.9877 

0.8687 

0.9206 

1.0007 

0.8332 

0.8076 

0.7956 

0.8236 

0.7999 

75 

VDOP 

2.3278 

2.3316 

2.3235 

2.3043 

1.8947 

1.9547 

2.2594 

1.745 

1.7581 

1.7713 

1.999 

2.0127 

VDOP/HDOP 

2.71 

2.72 

2.57 

2.33 

2.18 

2.12 

2.26 

2.09 

2.18 

2.23 

2.43 

2.52 
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